Avian lungs are remarkably different from mammalian lungs in that air flows unidirectionally through rigid tubes in which gas exchange occurs. Experimental observations have been able to determine the pattern of gas flow in the respiratory system, but understanding how the flow pattern is generated and determining the factors contributing to the observed dynamics remains elusive. It has been hypothesized that the unidirectional flow is due to aerodynamic valving during inspiration and expiration, resulting from the anatomical structure and the fluid dynamics involved, however, theoretical studies to back up this hypothesis are lacking. We have constructed a novel mathematical model of the airflow in the avian respiratory system that can produce unidirectional flow which is robust to changes in model parameters, breathing frequency and breathing amplitude. The model consists of two piecewise linear ordinary differential equations with lumped parameters and discontinuous, flow-dependent resistances that mimic the experimental observations. Using dynamical systems techniques and numerical analysis, we show that unidirectional flow can be produced by either effective inspiratory or effective expiratory valving, but that both inspiratory and expiratory valving are required to produce the high efficiencies of flows observed in avian lungs. We further show that the efficacy of the inspiratory and expiratory valving depends on airsac compliances and airflow resistances that may not be located in the immediate area of the valving. Our model provides additional novel insights; for example, we show that physiologically realistic resistance values lead to efficiencies that are close to maximum, and that when the relative lumped compliances of the caudal and cranial airsacs vary, it affects the timing of the airflow across the gas exchange area. These and other insights obtained by our study significantly enhance our understanding of the operation of the avian respiratory system.
The anatomical structure and airflow dynamics of the avian respiratory system are remarkably different to that of mammalian lungs [1]. The anatomical structure is complex, with multiple flexible airsacs that act like bellows to ventilate rigid tubes (parabronchi) in which gas exchange occurs, and a complicated branching structure that produces aerodynamic valving [1,2]. The airflow through the parabronchi (lungs) is unidirectional; flowing from the caudal (back) group of airsacs to the cranial (front) group of airsacs during both inspiration and expiration. (More precisely, the flow is unidirectional through the paleopulmonic-parabronchi that lie between the caudal and cranial airsacs. Some birds also contain neopulmonic-parabronchi in which airflow is bidirectional, but it forms a small part of the gas exchange surface area—less than 30% [2]. In this paper we use the term parabronchi to refer to the paleopulmonic parabronchi unless otherwise indicated.) Unlike in the mammalian respiratory system, the functions of ventilation and gas exchange have been uncoupled in the avian respiratory system; specifically, the flow of air through the system is caused by large flexible airsacs, whilst gas exchange occurs in narrow parabronchi which are rigid and firmly bound to the ribs [2]. The narrow, rigid structure of the parabronchi is thought to be related to the finding that birds have a thinner but mechanically stronger blood-gas barrier than equivalent mammals [3,4]. Furthermore the structure of the parabronchi and blood capillaries allows for cross-current gas exchange. These features are thought to contribute to the increased gas exchange efficiency of birds compared to mammals, especially at high-altitude or in a hypoxic environment [3,5–8]. The airflow pattern within the avian respiratory system is widely agreed upon. It has been determined by direct measurements of flow rates [9–11], as well as by experiments that used tracer gas, or CO2 and O2 measurements to indirectly determine the flow [10,12–15]. An important factor leading to unidirectional flow is hypothesized to be the effective inspiratory and expiratory aerodynamic valving that results from the interaction between the complex anatomical structure, including airway branching and constrictions, and the fluid dynamics involved [1,2]. The relative importance of the two valves has not been investigated. Additionally, the complex branching structure within the system affects the resistance to airflow of the different sections of the system. The effect of these resistances and the importance of their relative differences in generating the flow pattern is not known. Recently unidirectional airflow has been found in the lungs of some reptiles (specifically alligators [16,17], crocodiles [18,19], iguanas [20], and monitor lizards [21]). Comparing avian and reptile systems, which have very different levels of anatomical complexity in terms of the branching structures and the presence or absence of airsac separation, will provide important insights into the aerodynamic valving in both birds and reptiles. Mathematical modelling of the avian respiratory system has focused mainly on the gas exchange in the parabronchi (for example, [6,22–28]). These studies determined that gas exchange is cross-current and found gas exchange parameters for a range of avian species and experimental conditions [1,2]. Existing mathematical models of the airflow through the avian system had limited success in producing unidirectional flow [29,30]. Urushikubo et al. [30] used a three dimensional spatial model with simplified geometry for the pathways within the respiratory system, coupled with flexible airsacs. They found unidirectional flow through the parabronchi, but only for some parameter values. Additionally, the flow did not show any inspiratory or expiratory valving. Maina et al. [29] investigated aerodynamic inspiratory valving in ostriches by constructing a three dimensional anatomical model of the junction between the ventrobronchial branches and the main mesobronchus (i. e. the junction between the airways that lead air to the caudal airsacs and the airways that lead air from the cranial airsacs). Using computational fluid dynamics simulations they were only able to reproduce inspiratory valving if they included additional branches downstream (the secondary dorsobronchi branches), showing the importance of including the whole system when investigating aerodynamic valving. Unidirectional flow exists in all birds, despite massive inter-species differences in anatomy, and across most experimental conditions—including when ventilating the respiratory system post-mortem [31]. Thus, a useful mathematical model of the airflow in the avian respiratory system must produce unidirectional flow through the parabronchi across a broad range of parameter values and frequencies. In this paper we present a new, relatively simple, mathematical model of avian respiration that reproduces the airflow pattern described above. The unidirectional flow in our model is robust to changes in frequency and model parameters, and has efficiencies, flow rates, and pressures that match experimental findings. Additionally, our model generates several novel insights on the role of inspiratory and expiratory valving and the importance of variations in the airflow resistances and airsac compliances within the system that are thought to occur during respiration and in response to stimuli including hypoxia (lack of oxygen) and hypercapnia (excess of carbon dioxide). We first describe the mathematical model and then the new insights it produced. The model development and mathematical analysis are described later in the Methods section.
A schematic model of the avian respiratory system is shown in Fig 1. For simplicity, only one side of the respiratory system is shown (see Fig 12 for the full model). The caudal and cranial airsacs are considered to be flexible with lumped compliances C1 and C2 respectively and averaged pressures P1 and P2 respectively. The pressure in the coelom (thoracic-abdominal cavity) outside both sets of airsacs, Pext (t), varies periodically due to the respiratory muscles, which causes the airsacs to inflate and deflate. During inspiration (indicated by blue solid arrows in Fig 1), air flows in through the beak along the trachea (qT), through the primary and meso-bronchi to the caudal airsacs (q1), and from the caudal airsacs to the cranial airsacs through the parabronchi (qP). During expiration (indicated by green dashed arrows in Fig 1), air flows from the caudal airsacs to the cranial airsacs through the parabronchi (qP), from the cranial airsacs through the ventrobronchi (q2) and along the trachea to exit the beak (qT). The airflow pathways are considered to be rigid (no compliance) and to have resistance to airflow, Ri, where i ∈ {1,2, T, P}. Note that since our aim is to create a model that is applicable generally across all birds, and we are primarily interested in understanding the unidirectional flow through the paleopulmonic-parabronchi, we have chosen to include only the paleopulmonic-parabronchi explicitly. However, the neopulmonic-parabronchi are included indirectly in our model through the lumped resistance parameters. From the pressures P1, P2 and Patm (atmospheric pressure) and the resistances, Ri, we can calculate the airflow, qi, through every pathway in the system (Eqs (17) – (20) ). Using the relationship between the pressure, compliance, and volume, assuming that the compression of air is negligible, and applying some algebraic manipulations (refer to the Methods section for more details), we get the following two equations for the rate of change of the pressures P1 and P2: d P 1 d t = - R P R 2 (P 1 - P a t m) - (R P R T + R ¯) (P 1 - P 2) C 1 R ¯ R P + d P e x t d t (1) d P 2 d t = - R P R 1 (P 2 - P a t m) + (R P R T + R ¯) (P 1 - P 2) C 2 R ¯ R P + d P e x t d t (2) where R ¯ = R 1 R 2 + R 2 R T + R T R 1. The resistances R1 and R2 are discontinuous and vary depending on the flow direction as shown in Fig 2. This allows us to produce effective inspiratory and expiratory valving: A representative example of the results found in this model is shown in Fig 3, with the parameter values listed in Table 1. In Fig 3A we see that the pressure differences between the airsacs and atmospheric pressure (x1 = P1 − Patm and x2 = P2 − Patm) are orders of magnitude greater than the difference between the two airsac pressures (x1 − x2). This matches experimental measurements, e. g. [12,15]. Fig 3B shows that the flow through the parabronchi is unidirectional (qP > 0) and that the valving is working well as q2 ≈ 0 during inspiration and q1 ≈ 0 during expiration. The tidal volume is 36. 0 mL and the combined flow through the parabronchi per breath is 31. 3 mL on both sides, so most of the air which is breathed in passes through the gas exchange area. Fig 3C shows the volumes of the caudal set of airsacs, V1, and the cranial set of airsacs, V2, on one side of the respiratory system (see also Fig 12). The volumes can be calculated directly from Eqs (9) and (10). The ventilation volume into each set of airsacs is calculated using max (Vi) -min (Vi) for i = 1,2 and is independent of the chosen parameters Pc, and Vi, res. For the results shown in Fig 3 the ventilation volume is found to be 21. 5 mL per breath in total for the caudal airsacs, on both sides, and 16. 3 mL per breath in total for the cranial airsacs, on both sides. These values match experimental data [38]. Note that the sum of the ventilation of all the airsacs can be greater than or less than the tidal volume, as some air flows past the airsacs, and some air flows into both sets of airsacs. By analysing the phase plane dynamics of the system of Eq (1) (see Methods section), we can show that airflow through the parabronchi will be unidirectional (qP > 0) when γR1 ≤ R2 during inspiration and γR1 ≥ R2 during expiration, where γ = C1/C2. In the borderline case, where γR1 = R2 during both inspiration and expiration, the flow through the parabronchi, qP, is zero. A combination of effective inspiratory and expiratory valving (γR1 < R2 during inspiration and γR1 > R2 during expiration) will produce unidirectional flow. However, unidirectional flow could also be achieved by inspiratory or expiratory valving alone, e. g. effective inspiratory valving, where γR1 < R2 during inspiration and γR1 = R2 during expiration, or effective expiratory valving, where γR1 = R2 during inspiration and γR1 > R2 during expiration. Fig 4 shows a sketch of the phase plane dynamics in the case of effective inspiratory valving, where the variables are transformed such that x1 = P1 − Patm and x2 = P2 − Patm. The stable equilibrium P1 = P2 = Patm in the absence of pressure variations (no breathing) is then at the origin (0,0). The line x2 = x1 marks all the possible pressures for which the flow qP is zero. Above this line (x2 > x1) qP is negative (marked in shaded grey), and below it (x2 < x1) qP is positive. The dark blue curves show the solutions to the system from different initial conditions when there is no breathing, and thus no change in the external pressure. All these solutions approach the origin (0,0). The red curve shows the solution (not to scale) of the system when the external pressure Pext is changing due to breathing. This change in pressure is the same outside both sets of airsacs and thus acts along the vector [1,1]T. It can be seen that the system is ‘trapped’ below the line x1 = x2 and therefore the airflow is unidirectional (see Methods section for more details). In Fig 4 inspiration is marked by the light blue region and expiration is marked by the green region. When both P1 and P2 are greater than Patm (upper right quadrant) it is obvious that expiration will occur (qT < 0). When both P1 and P2 are less than Patm (lower left quadrant), inspiration will occur (qT > 0). The exact place (in the lower right quadrant) where the transition from inspiration to expiration occurs in the phase plane will depend on the chosen parameter values. Specifically, the flow qT will change direction on the line x 2 = - R 2, i n s p R 1, e x p x 1 with inspiration occurring if x 2 < - R 2, i n s p R 1, e x p x 1 and expiration occurring if x 2 > - R 2, i n s p R 1, e x p x 1 (see Methods section and Fig 13). For simplicity, in Fig 4 we consider the case where R2, insp = R1, exp and the transition between inspiration and expiration occurs on the line x2 = −x1. The conditions for unidirectional flow do not depend on the frequency or amplitude of breathing. Consequently, unidirectional flow will persist as long as the conditions on the resistances and compliances (stated in the previous section) are satisfied. Increasing the amplitude, Pamp, increases the flow rates proportionally (see Fig 5). As the period, T, decreases (the frequency increases) the mean flow through the parabronchi remains relatively constant, but the variation in the flow rate during the breathing cycle decreases and the flow becomes more constant (see Fig 6). This matches what is seen experimentally [11], and may have an impact on the efficacy of gas exchange. In the avian system, the aerodynamic valving is not 100% effective. Inspiratory valving is found experimentally to be 95–100% effective [14,32–34,39]. However, expiratory valving efficacy varies greatly between species and experimental conditions: 76–90% in chickens [12,40], 88% in ducks [14], and 95% in geese [37], with a strong dependence on gas velocity; at higher flow rates (exercise conditions) the valve is more effective than at rest. When fresh air flows into the cranial or caudal airsacs and is then breathed back out without passing through the parabronchi, this air does not undergo gas exchange, and is thus wasted. We define the efficiency of the whole system as the fraction of the tidal volume that passes through the parabronchi. For example, when efficiency = 1 all the air that is inhaled will pass through the parabronchi and undergo gas exchange. We calculate the efficiency in our model by numerically integrating the airflow qP during one cycle of breathing to find the total volume of air that flows through the parabronchi per breath, and numerically integrating the flow qT during inspiration to calculate the tidal volume (the total air inhaled per breath). The ratio of volume through parabronchi per breath to tidal volume gives us a measure of how efficient the lung system is. Efficiency = ∫ INSP + EXP q P dt ∫ INSP q T dt (5) where ∫INSP indicates the definite integral during inspiration and ∫INSP+EXP indicates the definite integral during one breath. If the model includes effective inspiratory valving only (R2, insp ≫ γR1, insp and R2, exp = γR1, exp) with R1, insp = R1, exp, the flow qP is unidirectional (qP > 0), but the maximum efficiency we can find numerically is around 50%. The cause of this low efficiency is that a large proportion of the fresh air that flows into the caudal airsacs, then flows back out (q1 < 0) without participating in gas exchange, as shown in Fig 7A. In Fig 7A, the inspiratory valving efficacy = 98. 7%, the expiratory valving efficacy = 47. 7%, the overall efficiency is 47. 1%, the tidal volume is 38. 1 mL whilst the flow through both parabronchi per breath is only 17. 9 mL. Similarly, if the model includes only effective expiratory valving (γR1, insp = R2, insp and γR1, exp ≫ R2, exp) with R2, insp = R2, exp, we find numerically that the maximum efficiency we can reach is around 50%, due to flow into the cranial airsacs during inspiration (q2 < 0). An example of effective expiratory valving is shown Fig 7B, where the inspiratory valving efficacy is 48. 6% and the expiratory valving efficacy is 87. 0%. This gives an overall efficiency of 42. 3%; from a tidal volume of 38. 4 mL and only 16. 2 mL flow through both parabronchi per breath. By including both inspiratory and expiratory valving we find that we can reduce this back flow, and when we match the experimental valving efficiencies of 98–100% for inspiratory valving and ≈88% for expiratory valving, we find that the overall efficiency of the system matches those found experimentally, whilst maintaining realistic resistance and compliance values. The flow rates for our chosen parameter values are shown in Fig 3B, where the inspiratory valving efficacy is 98. 0%, the expiratory valving efficacy is 88. 6%, and the overall efficiency is 86. 8%. In our model, we can investigate the impact of varying the resistances R1, insp and R2, exp. We need to keep R1, insp + R2, exp constant so that there is no change in the total resistance of the system, here we choose R1, insp + R2, exp = 6 cmH2O/L⋅s. Additionally, it is important to keep R2, insp = 100 × R1, insp and R1, exp = 10 × R2, exp, so that the strength of the valving isn’t changing. Fig 8 shows that depending on the ratio of compliances, γ, the maximum efficiency will be for 0. 2 < R1, insp/R2, exp < 1. This is consistent with experimental observations that found R1, insp to be lower than R2, exp (see Methods section). Comparatively, we find that varying γ and Ctot affect the efficiency by less than 1% (see Selecting model parameters). The airflow through the parabronchi is not constant during the breathing cycle. An important feature of the flow through the parabronchi is that it can be observed to occur mostly during inspiration, or expiration, or both, depending on parameter values and experimental conditions [10,11]. From Fig 7 we can see that the aerodynamic valving affects the timing of the airflow through the parabronchi; effective inspiratory valving increases parabronchial flow during inspiration (Fig 7A), whereas effective expiratory valving increases parabronchial flow during expiration (Fig 7B). However, once we fix the valving efficacy to physiologically realistic levels (e. g. inspiratory valving 98%, expiratory valving 88%), the ratio of volume flowing through the parabronchi during inspiration and expiration (I: E parabronchial volume flow ratio) due to the valving is fixed. We calculate the ratio of volume flowing through the parabronchi during inspiration and expiration from the output of our model as follows: I: E parabronchial volume flow ratio = ∫ INSP q P dt ∫ EXP q P dt (6) For the default parameter values (Table 1) the I: E parabronchial volume flow ratio is 0. 862, i. e. there is slightly less flow during inspiration than during expiration, as shown in Fig 3. Varying γ has a major impact on the timing of the flow through the parabronchi (recall that γ = C1/C2). When γ is low the majority of the flow through the parabronchi occurs during inspiration, while when γ is high the flow through the parabronchi occurs mostly during expiration. In Fig 9 we plot the I: E parabronchial volume flow ratio as a function of γ. As γ increases we find that the majority of the flow qP moves from being during the inspiratory phase to being during the expiratory phase. This result is conserved for a range of total compliance (C1 + C2 = Ctot) values. Despite this change in the timing of the flow, the system’s overall efficiency only decreases slightly (from 87. 7% to 86. 4%) when γ increases (see S1A Fig), and the airflow through the parabronchi remains unidirectional. These results are consistent with experimental observations. Anatomically, the caudal and cranial airsacs are found to have different properties [2]. In ducks, Scheid et al. [38] found that the caudal airsacs are more compliant and have larger ventilation volume changes than the cranial airsacs, especially during relaxed (anaesthetized) breathing. Furthermore, the ratio of compliances varies between individuals and species [1] and many variations in the flow pattern are also observed [10,11]. For example, in spontaneously breathing geese the flow through the parabronchi increases at the end of inspiration and peaks during expiration [10]. Looking at ducks it is found that the flow during inspiration is higher than the flow during expiration when panting, the flow during inspiration and expiration are similar in spontaneous breathing, and when relaxed (anaesthetized) the flow rate is much stronger during expiration [11]. Our results as we vary γ provides similar changes in flow patterns (Fig 10), which agree with experimental findings; γ should decrease during exercise when the abdominal and chest muscles are stiff, and increase under relaxation conditions when the muscles relax. We find that changing the relative resistances R1, insp/R2, exp does affect the timing of the flow through the parabronchi slightly, with more flow during expiration as R1, insp/R2, exp increases (see S2A Fig). We also find that varying the total compliance Ctot does not change the timing of the flow through the parabronchi substantially, but the strength of the effect decreases at high Ctot (see S2B Fig). Overall, the ratio of compliances, γ, is the dominant effect. We observe that although the forcing of the system is symmetric (sinusoidal function), the duration of inspiration, Ti (measured as the time during which qT > 0), is not always equal to the duration of expiration, Te (measured as the time during which qT < 0). For our chosen default parameter values (Table 1), with period T = 3s, Ti = 1. 4s and Te = 1. 6s, and the ratio of inspiration duration to expiration duration (I: E time ratio = Ti/Te) is 0. 89. This asymmetry varies depending on parameter values. We investigate the impact of varying the resistances R1, insp and R2, exp, while the overall resistance of the system and the strength of the valving constant, as before. We find that when we decrease R1, insp relative to R2, exp the duration of expiration increases, with a concordant decrease in the duration of inspiration (Fig 11). The ratio of compliances, γ, and the total compliance, Ctot, do affect the I: E time ratio slightly, but the impact of varying the relative resistances is much stronger (see Selecting model parameters).
The complex anatomical structure of the avian respiratory system has been represented in our model by discontinuous resistances (R1 and R2) that depend on the direction of airflow through them. These resistance values could depend on the properties of the flow and would then vary with frequency and amplitude of breathing, as well as with other parameters such as muscle tone that we did not take into account in our model. Nevertheless, we have shown that as long as γR1 ≤ R2 during inspiration and γR1 ≥ R2 during expiration, where γ = C1/C2 is the ratio of the airsac compliances, unidirectional flow will persist. We also assumed that both the inspiratory and expiratory valving are highly effective which is true during regular breathing but may not be the case if breathing consists of very high or very low frequencies or amplitudes. In particular, we note that experimentally it is found that panting and other breathing patterns (bird song/calls) do not have the same pattern. For example, during panting the expiratory valving is not strong and there is a large amount of air shunted into the primary bronchus (q1 < 0) which bypasses the parabronchi [9]. The two discontinuous resistances in our model make the system nonlinear, despite the assumptions of constant resistance and compliance elements. The presence of discontinuities in models is known to produce complicated phenomena, especially in non-autonomous systems with external forcing [41]. In this model, we have found the intriguing result that the inspiratory and expiratory periods are uneven in response to regular (sinusiodal) forcing. The phenomenon underlying this disparity is not clear yet. Further theoretical analysis is left for future investigations. In summary, the new mathematical model we have developed, and the analytical and computational study we have conducted, significantly increase our understanding of unidirectional airflow in avian lungs. Our new model is broadly applicable across all birds and can be extended or integrated into larger systems-level studies of the avian respiratory system. Our model also provides a new example of a non-smooth dynamical system and will be used in future investigations of the human respiratory system through comparative physiology.
Fig 12 shows the full model including left and right sides of the respiratory system. The caudal and cranial airsacs have pressures P1 and P2 respectively, and compliances C1 and C2 respectively. All other pathways and junctions are assumed to be rigid. To simplify our analysis, we assume that the left and right sides of the bird are symmetrical; this enables us to reduce the model to the system shown in Fig 1, where RT = 2Rtrachea + REPPB. In the remainder of this work we will use the model shown in Fig 1, which considers only one side. To find the overall flow in the whole animal, we simply double the flow rates found from this single side. To construct the mathematical model we begin by calculating the rate of change of volume in the caudal and cranial airsacs (dV1/dt and dV2/dt, respectively), assuming that the compression of air is negligible: d V 1 dt = q 1 - q P (7) d V 2 dt = q P - q 2 (8) where qi with i ∈ {1,2, P} is the airflow through the corresponding section. Next we assume that the airsacs are elastic, with compliance C1 and C2. Additionally, surrounding both sets of airsacs there is an external pressure Pext (t), which has a time-varying component that represents the change in pressure generated by the muscles of the chest and abdomen during breathing. This gives the equations: V 1 = C 1 P 1 - P e x t + V 1, r e s (9) V 2 = C 2 P 2 - P e x t + V 2, r e s (10) where V1, res and V2, res are the resting volumes of the caudal and cranial airsacs when the pressure difference between the airsacs and the surrounding thoracic-abdominal cavity (coelom) is zero. In all our simulations we use the sinusoidal function: P e x t (t) = P c - P a m p 2 cos 2 π t T (11) to model the time-varying pressure outside the airsacs. This function oscillates with a peak-to-peak amplitude of Pamp, which is the amplitude of the forcing from breathing, around a pressure Pc, which is the baseline pressure in the coelom. Differentiating Eqs (9) and (10) with respect to time gives d V 1 dt = C 1 d P 1 dt - d P e x t dt (12) d V 2 dt = C 2 d P 2 dt - d P e x t dt (13) with d P e x t dt = P a m p π T sin 2 π t T (14) Equating Eqs (12) and (7), and Eqs (13) and (8) gives: d P 1 dt = q 1 - q P C 1 + d P e x t dt (15) d P 2 dt = q P - q 2 C 2 + d P e x t dt (16) Assuming laminar flow, we obtain expressions for the flow rates qT, q1, q2, and qP in terms of our variables P1 and P2, as well as the pressures Patm and PJ: q T = P a t m - P J R T (17) q 1 = P J - P 1 R 1 (18) q 2 = P 2 - P J R 2 (19) q P = P 1 - P 2 R P (20) From the geometry of the system, the flow at junction J is conserved (inflexible junction), so qT + q2 = q1. From this, and flow Eqs (17) – (19) we find: R T (q 1 - q 2) + P J = P a t m (21) - R 1 q 1 + P J = P 1 (22) R 2 q 2 + P J = P 2 (23) Solving for q1, q2, and PJ in terms of P1, P2, and Patm we get: q 1 = - R 2 (P 1 - P a t m) - R T (P 1 - P 2) R ¯ (24) q 2 = R 1 (P 2 - P a t m) - R T (P 1 - P 2) R ¯ (25) P J = R 2 R T P 1 + R T R 1 P 2 + R 1 R 2 P a t m R ¯ (26) where we introduce the combined resistance R ¯ = R 1 R 2 + R 2 R T + R T R 1 to simplify the equations. Substituting the expressions for the flow rates Eqs (20), (24) and (25) into Eqs (15) and (16), we get the rate Eqs (1) and (2) for P1 and P2 respectively. The unique steady-state for the system of ordinary differential Eqs (1) and (2) in the absence of changing pressure due to breathing (d P e x t dt = 0), is P1 = P2 = Patm. If we move the equilibrium point to the origin using the transformation x1 = P1 − Patm, x2 = P2 − Patm, the system of equations becomes: d x 1 dt = - (R P R T + R P R 2 + R ¯) C 1 R ¯ R P x 1 + (R P R T + R ¯) C 1 R ¯ R P x 2 + d P e x t dt (27) d x 2 dt = (R P R T + R ¯) C 2 R ¯ R P x 1 + - (R P R T + R P R 1 + R ¯) C 2 R ¯ R P x 2 + d P e x t dt (28) The flow rates in terms of these new variables are: q T = - R 2 x 1 - R 1 x 2 R ¯ (29) q 1 = - R 2 x 1 - R T (x 1 - x 2) R ¯ (30) q 2 = R 1 x 2 - R T (x 1 - x 2) R ¯ (31) q P = x 1 - x 2 R P (32) Eqs (27) and (28) can be written in matrix form as: d X → dt = A X → + d P e x t dt 1 1 (33) where A = - (R P R T + R P R 2 + R ¯) C 1 R ¯ R P (R P R T + R ¯) C 1 R ¯ R P (R P R T + R ¯) C 2 R ¯ R P - (R P R T + R P R 1 + R ¯) C 2 R ¯ R P (34) and X → = [ x 1, x 2 ] T. If we consider the unforced system, where d P e x t dt = 0, Eq (33) becomes the autonomous, homogeneous, linear system d X U → dt = A X U →, where X U → = [ x 1, x 2 ] T. The solution to this autonomous linear system is: X U → (t) = a 1 e λ 1 t V 1 → + a 2 e λ 2 t V 2 → (35) where λi are the eigenvalues of A, V i → are their corresponding eigenvectors, and ai depend on the specific initial conditions. To find explicit expressions for the eigenvalues and eigenvectors of the unforced system, we simplify our analysis by scaling time by C 1 R ¯ R P. This transformation does not change the phase plane dynamics of the system. The matrix A (Eq (34) ) becomes: A ^ = - (R 2 R P + β) β γ β - γ (R 1 R P + β) (36) where β = R T R P + R ¯, and γ = C1/C2. The eigenvalues are then given by: λ 1 = T r (A ^) - T r (A ^) 2 - 4 Det (A ^) 2 (37) λ 2 = T r (A ^) + T r (A ^) 2 - 4 Det (A ^) 2 (38) where: T r (A ^) = - (γ R 1 + R 2) R P - (γ + 1) β (39) Det (A ^) = γ (R 1 R 2 R P 2 + β (R 1 + R 2) R P) (40) These eigenvalues will be real if Δ = Tr (A ^) 2 - 4 Det (A ^) ⩾ 0. In this system all resistance and compliance values must be positive. Using this simple physiological constraint we find that T r (A ^) < 0 and Det (A ^) > 0 for all values of Ci, Ri > 0. Additionally, it can be shown that Δ = Tr (A ^) 2 - 4 Det (A ^) > 0, and thus the system has two real eigenvalues, λ1 < λ2 < 0, for all parameter values. The corresponding eigenvectors are: V 1 → = 1 λ 1 + R P R 2 + β β (fast) (41) V 2 → = 1 λ 2 + R P R 2 + β β (slow) (42) All solutions (given by Eq (35) ) will tend to the single steady-state (x1, x2) = (0,0) as t → ∞ by initially following the fast eigenvector V 1 → and then approaching the steady-state (more slowly) following the slow eigenvector V 2 →. This behaviour is characteristic of solutions to linear ordinary differential equations. In our model, we are looking for solutions where qP ≥ 0. From Eq (32) we know that qP = 0 on the line x2 = x1 and solutions above the line x2 = x1 will have qP < 0 while solutions below the line x2 = x1 will have qP > 0. Thus, to maintain positive unidirectional flow (qP ≥ 0), solutions must lie on or below the line x2 = x1. Putting this information together with the knowledge that the unforced system will return to steady state (x1, x2) = (0,0) following the slow eigenvector V 2 →, we can conclude that solutions will be pushed into the region qP < 0 if the slow eigenvector, V 2 → lies above the line x2 = x1. Consequently, in order to maintain unidirectional flow we must find conditions that ensure that the slow eigenvector lies on or below the line x2 = x1. From Eq (42), V 2 → = [ 1,1 ] T when λ 2 + R P R 2 + β β = 1. Rearranging this equation, we find that the slow eigenvector will lie along the line x2 = x1 when γR1 = R2. When perturbed by breathing, Pext is applied along the vector [1,1]T (see Eq (33) ). So in the case where γR1 = R2 and V 2 → = [ 1,1 ] T, the system will oscillate (due to the forcing from breathing) along the line x2 = x1 and qP = 0. This forms a useful boundary case. During inspiration the slow eigenvector will lie below the line x1 = x2 (qP > 0) if λ 2 + R P R 2 + β β > 1, which occurs when γR1 < R2. During expiration, the slow eigenvector will lie below the line x1 = x2 (qP > 0) if λ 2 + R P R 2 + β β < 1, which occurs when γR1 > R2. If both these conditions are satisfied, then the dynamics of the unforced system tells us that solutions will move quickly into the region qP > 0 and will stay there. When the system is perturbed by breathing, the forcing Pext is applied along the vector [1,1]T, and thus will not move the solutions into the region qP < 0. This analysis gives the following conditions for unidirectional flow: γR1 ≤ R2 during inspiration and γR1 ≥ R2 during expiration. The conditions for unidirectional flow stated above can only be satisfied if the resistances or compliances change between inspiration and expiration. As there is no evidence that the compliances would change during a single breathing cycle, we look instead at changing the resistances. Experimental findings have shown that the complicated anatomical structure causes effective aerodynamic valving, where the flow through the different pathways differs between inspiration and expiration. To reproduce this effect in our model we make the resistances R1 and R2 dependent on flow direction. Inspiratory valving is incorporated into our model by increasing R2 when q2 < 0 (see Fig 2) to reduce the negative q2 flow during inspiration. The resistance R2 is: R 2 = R 2, exp + (R 2, i n s p − R 2, exp) H (P J − P 2) where H denotes the Heaviside function, R2, exp is the physiological value for resistance to flow in the preferred direction (q2 > 0) and R2, insp is the higher effective resistance value during inspiration due to the inspiratory valving. Expiratory valving is incorporated into our model by increasing R1 when q1 < 0 to prevent flow from the caudal airsacs during expiration (see Fig 2). The resistance R1 is: R 1 = R 1, i n s p + (R 1, exp − R 1, i n s p) H (P 1 − P J) where H denotes the Heaviside function, R1, insp is the physiological value for resistance to flow in the preferred direction (q1 > 0) and R1, exp is the higher effective resistance value during inspiration due to the expiratory valving. From Eq (30) the flow q1 = 0 when x 2 = R 2 + R T R T x 1, for qP > 0 this line lies in the lower left quadrant of the (x1, x2) -phase plane. From Eq (31) the flow q2 = 0 when x 2 = R T R 1 + R T x 1, for qP > 0 this line lies in the upper right quadrant of the (x1, x2) -phase plane. From Eq (29) we can show that inspiration (qT > 0) occurs in the region x 2 < - R 2 R 1 x 1, and expiration (qT < 0) occurs in the region x 2 > - R 2 R 1 x 1. Combining all this information we can sketch the flow direction transitions onto the phase plane (Fig 13), where: Note that the qT = 0 transition always occurs in the lower right quadrant, where q1 < 0 and q2 < 0. This means that we can define inspiration as the region x 2 < - R 2, i n s p R 1, exp x 1, and expiration as the region x 2 > - R 2, i n s p R 1, exp x 1. Also note that the transition through regions (2) and (3) is very fast. We used Matlab’s event detection in conjunction with the solver ode23 to change R1 and R2 values when q1 and q2 crossed through zero. Starting at region (1) (q1 > 0, q2 < 0) the sequence for one breath is: This sequence was then repeated for as many breaths as required until steady state was reached. Steady state was numerically defined as being reached when the area under the curve qT was zero (less than 1 × 10−5) over a single breath. That is, the flow into the bird was equal to the flow out of the bird in each breath. The volumetric flow (area under the qi curves) through each segment was found numerically using the trapezoidal method. In the results presented, a step size of δt = 1 × 10−4 was used. The step size was reduced in several cases to test convergence and the results were consistent. Note: We state and discuss resistances with the units cmH2O/L⋅s, but use the units mL/cmH2O for compliances and mL for volumes, based on common practice in the field. When implementing the model it is important to use consistent units (mL or L only). We selected parameters to match duck respiratory systems, as we have the best data on airsac compliance and ventilation for this species. The default parameter values given in Table 1 are used for all the numerical calculations unless otherwise is indicated in the figure legends or in the text. Below we explain in more detail how we chose the specific parameters.
Birds and mammals have similar metabolic demands and cardiovascular systems, but they have evolved drastically different respiratory systems. A key difference in birds is that gas exchange occurs in rigid tubes, through which air flows unidirectionally during both inspiration and expiration. How this unidirectional flow is generated, and the factors affecting it, are not well understood. It has been hypothesized that the unidirectional flow is due to aerodynamic valving resulting from the complex anatomical structure. To test this hypothesis we have constructed a novel mathematical model that, unlike previous models, produces unidirectional flow through the lungs consistently even when the amplitude and frequency of breathing change. We have investigated the model both analytically and computationally and shown the importance of aerodynamic valving for generating strong airflow through the lungs. Our model also predicts that the timing of airflow through the lungs depends on the relative compliances of the different airsacs that exist in birds. The lumped parameters approach we use means that this model is generally applicable across all birds.
livestock
medicine and health sciences
dynamics
classical mechanics
fluid mechanics
vertebrates
animals
physiological processes
respiratory system
breathing
respiration
waterfowl
air flow
birds
ducks
fluid dynamics
fowl
continuum mechanics
poultry
agriculture
fluid flow
physics
flow rate
anatomy
physiology
biology and life sciences
physical sciences
organisms
aerodynamics
Robust Unidirectional Airflow through Avian Lungs: New Insights from a Piecewise Linear Mathematical Model