Robust Unidirectional Airflow through Avian Lungs: New Insights from a Piecewise Linear Mathematical Model

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.


Introduction
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][6][7][8].
The airflow pattern within the avian respiratory system is widely agreed upon. It has been determined by direct measurements of flow rates [9][10][11], as well as by experiments that used tracer gas, or CO 2 and O 2 measurements to indirectly determine the flow [10,[12][13][14][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][23][24][25][26][27][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.

The model
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 C 1 and C 2 respectively and averaged pressures P 1 and P 2 respectively. The pressure in the coelom (thoracic-abdominal cavity) outside both sets of airsacs, P ext (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 (q T ), through the primary and mesobronchi to the caudal airsacs (q 1 ), and from the caudal airsacs to the cranial airsacs through the parabronchi (q P ). During expiration (indicated by green dashed arrows in Fig 1), air flows from the caudal airsacs to the cranial airsacs through the parabronchi (q P ), from the cranial airsacs through the ventrobronchi (q 2 ) and along the trachea to exit the beak (q T ). The airflow pathways are considered to be rigid (no compliance) and to have resistance to airflow, R i , where i 2 {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 paleopulmonicparabronchi explicitly. However, the neopulmonic-parabronchi are included indirectly in our model through the lumped resistance parameters.
From the pressures P 1 , P 2 and P atm (atmospheric pressure) and the resistances, R i , we can calculate the airflow, q i , 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 P 1 and P 2 : where The resistances R 1 and R 2 are discontinuous and vary depending on the flow direction as shown in Fig 2. This allows us to produce effective inspiratory and expiratory valving: • Inspiratory valving. During inspiration, it is observed that air flows through the junction P J into the caudal airsacs (q 1 % q T ), with very little fresh air flowing into the cranial airsacs (q 2 ( q T ). This is attributed to anatomical features including the T-shape of the junction, the narrowing of the airway at the segmentum accelerans, and the branching nature of the airway tree [2,14,29,[32][33][34][35][36]. We incorporate this valving effect into our model by increasing R 2 during inspiration (see Fig 2B) and we measure the effectiveness of the valving by calculating how much of the flow into the animal flows through to the caudal airsacs during inspiration Schematic model of the avian respiratory system. The caudal and cranial airsacs have pressures P 1 and P 2 respectively, and compliances C 1 and C 2 respectively. The pressure outside both sets of airsacs, P ext (t), varies periodically due to the respiratory muscles, which causes the airsacs to inflate and deflate. The pressure P atm is atmospheric pressure. Between each node in the system there is resistance to flow, R i , and airflow, q i . Blue arrows represent the flow during inspiration, and green arrows represent the flow during expiration. The grey shaded area indicates the parabronchi, where gas exchange occurs.
• Expiratory valving. During expiration, it is observed that most of the air from the caudal airsacs flows through the parabronci into the cranial airsacs (q P ), with very little fresh air flowing back without undergoing gas exchange (q 1 ( q T ). Expiratory valving is thought to be due to the specific anatomical structure and alignment of airsacs and dorsobronchial airways [37] but is not as well studied as the inspiratory valving. We incorporate this valving into our model by increasing R 1 during expiration (see Fig 2A) and we define the effectiveness of the valving as how much of the flow out of the system comes from the cranial airsacs during expiration (labeled EXP).
Expiratory valving efficacy ¼ Note that other definitions of expiratory valving efficacy exist in the literature, see Eq (49).

Model outputs
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 The resistance R 1 is discontinuous at P 1 = P J , which is when The resistance R 2 is discontinuous at P 2 = P J , which is when and atmospheric pressure (x 1 = P 1 − P atm and x 2 = P 2 − P atm ) are orders of magnitude greater than the difference between the two airsac pressures (x 1 − x 2 ). This matches experimental measurements, e.g. [12,15]. Fig 3B shows that the flow through the parabronchi is unidirectional (q P > 0) and that the valving is working well as q 2 % 0 during inspiration and q 1 % 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, V 1 , and the cranial set of airsacs, V 2 , 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 (V i )-min(V i ) for i = 1, 2 and is independent of the chosen parameters P c , and V i,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.  (Table 1). A: the pressure differences from atmospheric pressure, x 1 = P 1 − P atm and x 2 = P 2 − P atm , which are the outputs of the model. B: the flow rates q T , q P , q 1 , and q 2 in one side of the respiratory system. C: the volumes of the caudal set of airsacs, V 1 , and the cranial set of airsacs, V 2 , on one side of the respiratory system. Inspiration and expiration phases are labelled (INSP and EXP).

Conditions for unidirectional airflow
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 (q P > 0) when γR 1 R 2 during inspiration and γR 1 ! R 2 during expiration, where γ = C 1 /C 2 . In the borderline case, where γR 1 = R 2 during both inspiration and expiration, the flow through the parabronchi, q P , is zero. A combination of effective inspiratory and expiratory valving (γR 1 < R 2 during inspiration and γR 1 > R 2 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 γR 1 < R 2 during inspiration and γR 1 = R 2 during expiration, or effective expiratory valving, where γR 1 = R 2 during inspiration and γR 1 > R 2 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 x 1 = P 1 − P atm and x 2 = P 2 − P atm . The stable equilibrium P 1 = P 2 = P atm in the absence of pressure variations (no breathing) is then at the origin (0, 0). The line x 2 = x 1 marks all the possible pressures for which the flow q P is zero. Above this line (x 2 > x 1 ) q P is negative (marked in shaded grey), and below it (x 2 < x 1 ) q P 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 P ext 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 x 1 = x 2 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 P 1 and P 2 are greater than P atm (upper right quadrant) it is obvious that expiration will occur (q T < 0). When both P 1 and P 2 are less than P atm (lower left quadrant), inspiration will occur (q T > 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 q T will change direction on the line x 2 ¼ À

Unidirectional flow is robust to changes in the breathing amplitude and frequency
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, P amp , 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  . This matches what is seen experimentally [11], and may have an impact on the efficacy of gas exchange.
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 q P during one cycle of breathing to find the total volume of air that flows through the parabronchi per breath, and numerically integrating the flow q T 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  Table 1.  Unidirectional flow is robust to changes in frequency. The flow through the parabronchi, q P , during a single breath, is plotted for breathing period, T, ranging from 1-6 seconds (frequency ranging from 1-1/6 Hz respectively). The traces are aligned such that phase = 0 is at the beginning of inspiration. Inspiration and expiration are labelled (INSP and EXP). All the parameters, except the breathing period, are as in Table 1.
where R INSP indicates the definite integral during inspiration and R INSP+EXP indicates the definite integral during one breath.
If the model includes effective inspiratory valving only (R 2,insp ) γR 1,insp and R 2,exp = γR 1, exp ) with R 1,insp = R 1,exp , the flow q P is unidirectional (q P > 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 (q 1 < 0) without participating in gas exchange, as shown in Fig 7A. In Fig 7A,  Similarly, if the model includes only effective expiratory valving (γR 1,insp = R 2,insp and γR 1,exp ) R 2,exp ) with R 2,insp = R 2,exp , we find numerically that the maximum efficiency we can reach is around 50%, due to flow into the cranial airsacs during inspiration (q 2 < 0). An example of Inspiratory and expiratory valving both produce unidirectional flow, q P > 0. Flow rates q T , q P , q 1 , and q 2 against time for the parameters R 1,insp = R 2,exp = 3 cmH 2 O/LÁs, C 1 = C 2 (γ = 1), and C tot = 450 mL/cmH 2 O. Panel A shows the case where there is effective inspiratory valving: R 2,exp = γR 1,exp and R 2, insp = 100 × R 1,insp , with R 1,exp = R 1,insp . Panel B shows the case with effective expiratory valving: R 1,insp = γR 2,insp and R 1,exp = 20 × R 2,exp , with R 2,insp = R 2,exp . All other parameters are given in Table 1. 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%.
Efficiency is affected more by asymmetries in the resistance to flow than by the compliances In our model, we can investigate the impact of varying the resistances R 1,insp and R 2,exp . We need to keep R 1,insp + R 2,exp constant so that there is no change in the total resistance of the system, here we choose R 1,insp + R 2,exp = 6 cmH 2 O/LÁs. Additionally, it is important to keep R 2,insp = 100 × R 1,insp and R 1,exp = 10 × R 2,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 < R 1,insp /R 2,exp < 1. This is consistent with experimental observations that found R 1,insp to be lower than R 2,exp (see Methods section). Comparatively, we find that varying γ and C tot affect the efficiency by less than 1% (see Selecting model parameters).
The timing of the airflow through the parabronchi depends on the relative airsac compliance, γ 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]. Changing the relative resistance of R 1,insp / R 2,exp affects the efficiency of the system. Plot of the overall efficiency when R 1,insp / R 2,exp is varied whilst keeping the total resistance of the system constant (R 1, insp + R 2,exp = 6 cmH 2 O/LÁs). The effect is similar for a range of γ values. All other parameters are given in Table 1.
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: 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 γ = C 1 /C 2 ). 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 q P moves from being during the inspiratory phase to being during the expiratory phase. This result is conserved for a range of total compliance (C 1 + C 2 = C tot ) 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 R 1,insp /R 2,exp does affect the timing of the flow through the parabronchi slightly, with more flow during expiration as R 1,insp /R 2,exp increases (see S2A Fig). We also find that varying the total compliance C tot does not change the timing of the flow through the parabronchi substantially, but the strength of the effect decreases at high C tot (see S2B Fig). Overall, the ratio of compliances, γ, is the dominant effect.
The durations of inspiration and expiration depend primarily on the relative resistances R 1,insp and R 2,exp We observe that although the forcing of the system is symmetric (sinusoidal function), the duration of inspiration, T i (measured as the time during which q T > 0), is not always equal to the duration of expiration, T e (measured as the time during which q T < 0). For our chosen default parameter values (Table 1), with period T = 3s, T i = 1.4s and T e = 1.6s, and the ratio of inspiration duration to expiration duration (I:E time ratio = T i /T e ) is 0.89. This asymmetry varies depending on parameter values.
We investigate the impact of varying the resistances R 1,insp and R 2,exp , while the overall resistance of the system and the strength of the valving constant, as before. We find that when we decrease R 1,insp relative to R 2,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, C tot , do affect the I:E time ratio slightly, but the impact of varying the relative resistances is much stronger (see Selecting model parameters).

Discussion
We have constructed a relatively simple mathematical model of the avian respiratory system which, for the first time, produces unidirectional airflow through the parabronchi that is robust to changes in breathing frequency, breathing amplitude, and model parameters. Using   Fig 10. The ratio of compliances γ = C 1 /C 2 changes the shape of the oscillatory flow q P . This figure plots the flow rate q P versus time, for γ = 1/4 (red), γ = 1 (green), and γ = 4 (blue). The inspiratory period (INSP) is shaded blue, and the expiratory period (EXP) is shaded green. All other parameters are as given in Table 1. physiologically reasonable parameters (in this case we use those found for ducks) the model produces efficiencies, flow rates, and pressures that match experimental findings.
It has been hypothesized that the unidirectional flow is due to inspiratory and expiratory aerodynamic valves, resulting from the anatomical structure and the fluid dynamics involved. We incorporated these aerodynamic valves into our model by increasing the resistance to flow from the primary bronchus to the cranial airsacs (R 2 in our model) during inspiration, and increasing the resistance to flow from the caudal airsacs to the primary bronchus (R 1 in our model) during expiration. We showed that both of these resistances as well as the airsac compliances (C 1 and C 2 in our model) affect the efficacy of the inspiratory and expiratory valving. This result could explain why models that focused on limited areas of the respiratory system or that oversimplified the pathways' geometry were unable to produce the valving [29,30]. We further showed that unidirectional flow could be produced by either an effective inspiratory or an effective expiratory valve (Fig 7), but that both inspiratory and expiratory valves are required to produce the high efficiencies observed in avian lungs.
In existing models, the compliances of the caudal and cranial airsacs has been assumed to be the same [29,30]. However, there is no anatomical reason why this should be the case, and indeed this is not what has been found experimentally [38]. Using our model, we varied the single parameter γ = C 1 /C 2 , whilst keeping the total compliance constant (C tot = C 1 + C 2 ). We found that the ratio of compliances does not affect the total flow through the parabronchi significantly, but that it has a strong impact on the timing of the flow; for physiological parameter values when C 1 < C 2 , the majority of the flow through the parabronchi occurs during inspiration, and when C 1 > C 2 , the majority of the flow through the parabronchi occurs during expiration (Fig 9). The overall compliance of the airsacs is affected by the chest wall and muscles surrounding them. This means that the effective compliance is a parameter that could vary in different conditions, and would strongly influence the dynamics of the flow through the parabronchi and thus the gas exchange.
Our model provides additional novel insights into the operation of the avian respiratory system. We showed that changing the relative resistance of R 1,insp / R 2,exp whilst keeping the total resistance (R 1,insp + R 2,exp ) constant, affects the efficiency of the system and that maximum efficiencies appear to exist near physiologically realistic parameter values (Fig 8). Another interesting observation is that the ratio of R 1,insp /R 2,exp affects the period of expiration and inspiration. Specifically, we found that T e > T i when R 1,insp < R 2,exp (Fig 11).

Limitations
The complex anatomical structure of the avian respiratory system has been represented in our model by discontinuous resistances (R 1 and R 2 ) 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 γR 1 R 2 during inspiration and γR 1 ! R 2 during expiration, where γ = C 1 /C 2 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 (q 1 < 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.

Conclusions
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 P 1 and P 2 respectively, and compliances C 1 and C 2 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 R T = 2R trachea + R EPPB . 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.

Model formulation
To construct the mathematical model we begin by calculating the rate of change of volume in the caudal and cranial airsacs (dV 1 /dt and dV 2 /dt, respectively), assuming that the compression of air is negligible: where q i with i 2 {1, 2, P} is the airflow through the corresponding section. Next we assume that the airsacs are elastic, with compliance C 1 and C 2 . Additionally, surrounding both sets of airsacs there is an external pressure P ext (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: where V 1,res and V 2,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: to model the time-varying pressure outside the airsacs. This function oscillates with a peak-topeak amplitude of P amp , which is the amplitude of the forcing from breathing, around a pressure P c , which is the baseline pressure in the coelom. Differentiating Eqs (9) and (10) with respect to time gives Equating Eqs (12) and (7), and Eqs (13) and (8) gives: Assuming laminar flow, we obtain expressions for the flow rates q T , q 1 , q 2 , and q P in terms of our variables P 1 and P 2 , as well as the pressures P atm and P J : From the geometry of the system, the flow at junction J is conserved (inflexible junction), so q T + q 2 = q 1 . From this, and flow Eqs (17)- (19) we find: Solving for q 1 , q 2 , and P J in terms of P 1 , P 2 , and P atm we get: where we introduce the combined resistance 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 P 1 and P 2 respectively.

Analysis of the model
The unique steady-state for the system of ordinary differential Eqs (1) and (2) in the absence of changing pressure due to breathing dP ext dt ¼ 0 À Á , is P 1 = P 2 = P atm . If we move the equilibrium point to the origin using the transformation x 1 = P 1 − P atm , x 2 = P 2 − P atm , the system of equations becomes: The flow rates in terms of these new variables are: Eqs (27) and (28) can be written in matrix form as: If we consider the unforced system, where dP ext dt ¼ 0, Eq (33) becomes the autonomous, homogeneous, linear system dX U dt ¼ AX U , whereX U ¼ ½x 1 ; x 2 T . The solution to this autonomous linear system is:X where λ i are the eigenvalues of A,Ṽ i are their corresponding eigenvectors, and a i 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 RR P . This transformation does not change the phase plane dynamics of the system. The matrix A (Eq (34)) becomes: where b ¼ R T R P þ R, and γ = C 1 /C 2 . The eigenvalues are then given by: where: These eigenvalues will be real if D ¼ TrðÂÞ 2 À 4DetðÂÞ ⩾ 0.
In this system all resistance and compliance values must be positive. Using this simple physiological constraint we find that TrðÂÞ < 0 and DetðÂÞ > 0 for all values of C i , R i > 0. Additionally, it can be shown that D ¼ TrðÂÞ 2 À 4DetðÂÞ > 0, and thus the system has two real eigenvalues, λ 1 < λ 2 < 0, for all parameter values.
The corresponding eigenvectors are: All solutions (given by Eq (35)) will tend to the single steady-state (x 1 , x 2 ) = (0, 0) as t ! 1 by initially following the fast eigenvectorṼ 1 and then approaching the steady-state (more slowly) following the slow eigenvectorṼ 2 . This behaviour is characteristic of solutions to linear ordinary differential equations.
In our model, we are looking for solutions where q P ! 0. From Eq (32) we know that q P = 0 on the line x 2 = x 1 and solutions above the line x 2 = x 1 will have q P < 0 while solutions below the line x 2 = x 1 will have q P > 0. Thus, to maintain positive unidirectional flow (q P ! 0), solutions must lie on or below the line x 2 = x 1 . Putting this information together with the knowledge that the unforced system will return to steady state (x 1 , x 2 ) = (0, 0) following the slow eigenvec-torṼ 2 , we can conclude that solutions will be pushed into the region q P < 0 if the slow eigenvector,Ṽ 2 lies above the line x 2 = x 1 . Consequently, in order to maintain unidirectional flow we must find conditions that ensure that the slow eigenvector lies on or below the line Rearranging this equation, we find that the slow eigenvector will lie along the line x 2 = x 1 when γR 1 = R 2 . When perturbed by breathing, P ext is applied along the vector [1, 1] T (see Eq (33)). So in the case where γR 1 = R 2 and V 2 ¼ ½1; 1 T , the system will oscillate (due to the forcing from breathing) along the line x 2 = x 1 and q P = 0. This forms a useful boundary case. During inspiration the slow eigenvector will lie below the line x 1 = x 2 (q P > 0) if l 2 þR P R 2 þb b > 1, which occurs when γR 1 < R 2 . During expiration, the slow eigenvector will lie below the line x 1 = x 2 (q P > 0) if l 2 þR P R 2 þb b < 1, which occurs when γR 1 > R 2 . If both these conditions are satisfied, then the dynamics of the unforced system tells us that solutions will move quickly into the region q P > 0 and will stay there. When the system is perturbed by breathing, the forcing P ext is applied along the vector [1, 1] T , and thus will not move the solutions into the region q P < 0. This analysis gives the following conditions for unidirectional flow: γR 1 R 2 during inspiration and γR 1 ! R 2 during expiration.
Special case of even compliances. In the special case γ = 1 (that is, C 1 = C 2 ), the eigenvalues reduce to: and the corresponding eigenvectors are: If we set R 1 = R 2 , the eigenvalues simplify further to: and the eigenvectors are:Ṽ When perturbed by breathing, this system will return along the line x 1 = x 2 , where q P = 0, which forms the boundary (zero flow) case. If R 2 > R 1 , the slow eigenvectorṼ 2 (Eq (45)) lies below the line x 1 = x 2 during inspiration and q P > 0. If R 1 > R 2 ,Ṽ 2 (Eq (45)) lies below the line x 1 = x 2 during expiration and once again q P > 0. Fig 4 shows a representative phase plane for γ = 1.

Incorporating aerodynamic valving
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 R 1 and R 2 dependent on flow direction.
Inspiratory valving is incorporated into our model by increasing R 2 when q 2 < 0 (see Fig 2) to reduce the negative q 2 flow during inspiration. The resistance R 2 is: where H denotes the Heaviside function, R 2,exp is the physiological value for resistance to flow in the preferred direction (q 2 > 0) and R 2,insp is the higher effective resistance value during inspiration due to the inspiratory valving.
Expiratory valving is incorporated into our model by increasing R 1 when q 1 < 0 to prevent flow from the caudal airsacs during expiration (see Fig 2). The resistance R 1 is: where H denotes the Heaviside function, R 1,insp is the physiological value for resistance to flow in the preferred direction (q 1 > 0) and R 1,exp is the higher effective resistance value during inspiration due to the expiratory valving.
From Eq (30) the flow q 1 = 0 when x 2 ¼ R 2 þR T R T x 1 , for q P > 0 this line lies in the lower left quadrant of the (x 1 , x 2 )-phase plane. From Eq (31) the flow q 2 = 0 when x 2 ¼ R T R 1 þR T x 1 , for q P > 0 this line lies in the upper right quadrant of the (x 1 , x 2 )-phase plane. From Eq (29) we can show that inspiration (q T > 0) occurs in the region x 2 < À R 2 R 1 x 1 , and expiration (q T < 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: • In region (1): q T > 0, q 1 > 0, and q 2 < 0, so R 1 = R 1,insp , R 2 = R 2,insp .
Note that the q T = 0 transition always occurs in the lower right quadrant, where q 1 < 0 and q 2 < 0. This means that we can define inspiration as the region x 2 < À R 2;insp R 1;exp x 1 , and expiration as the region x 2 > À R 2;insp R 1;exp x 1 . Also note that the transition through regions (2) and (3) is very fast.

Implementation
We used Matlab's event detection in conjunction with the solver ode23 to change R 1 and R 2 values when q 1 and q 2 crossed through zero. Starting at region (1) (q 1 > 0, q 2 < 0) the sequence for one breath is: • set R 1 = R 1,insp , and R 2 = R 2,insp , and solve until q 1 = 0 (region 1), • set R 1 = R 1,exp , leave R 2 = R 2,insp , and solve until q 2 = 0 (regions 2 & 3), • leave R 1 = R 1,exp , set R 2 = R 2,exp , and solve until q 2 = 0 (region 4), • leave R 1 = R 1,exp , set R 2 = R 2,insp , and solve until q 1 = 0 (regions 2 & 3). 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 q T 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 q i 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 cmH 2 O/LÁs, but use the units mL/ cmH 2 O 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).

Selecting model parameters
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.
Resistances. R 2,exp = 5 cmH 2 O/LÁs and R P = 2.5 cmH 2 O/LÁs are chosen to match physiological values measured in ducks [42]. There are no direct measurements of R 1,insp , but we select R 1,insp = 1 cmH 2 O/LÁs to match measurements of the overall lower respiratory system resistance in a single side (R 1,insp + R P + R 2,exp % 9 cmH 2 O/LÁs) [42,43] and measurements that show that the pressure drop between the primary bronchi and the caudal airsacs is very small [34]. The tracheal resistance R trachea = 1 cmH 2 O/LÁs is chosen to match the resistance obtained from measurements in [43] and checked against the resistance to laminar flow in a rigid tube calculated for the approximate length and diameter [1,2,44,45]. The resistance across the constriction known as the segmentum accelerans R EPPB = 8 cmH 2 O/LÁs is chosen to generate observed pressure drop [34]. The overall respiratory system resistance can be calculated as R total = 2R trachea + R EPPB + (R 1,insp + R P + R 2,exp )/2. For our chosen parameter values R total = 14.25 cmH 2 O/LÁs, which is consistent with measurements [42,43,46,47].
Valving. To implement valving, we chose R 2,insp = 100 × R 1,insp to produce inspiratory valving efficacy in physiological range %96 − 100% [14,33,34], and chose R 1,exp = 10 × R 2,exp to produce expiratory valving efficacy in physiological range %88% for ducks at rest [14,48]. We note here that measurements of expiratory valving efficacy are often calculated using the amount of air that flows through the lungs q P as a proportion of the total amount of air that flows out of the caudal airsacs during expiration: This differs slightly from our definition of expiratory valving efficacy (see Eq (4)), and will give a lower estimate of the valving efficacy compared to our definition. However, using the above definition, our chosen parameters give 81% expiratory valving efficacy, which is still in the experimentally measured range of 76-95%. Total airsac compliance. Because of experimental limitations, the compliance of the airsacs has not been measured; measurements have only been made of the total system, including the body wall. One technique that is used for measuring the compliance is to change the external pressure in a box around an anaesthetized bird and measure the resulting volume flow into or out of the animal. The slope ΔV/ΔP then gives a measure of the 'static compliance' of the total system. Experiments in chickens and ducks estimate the compliance at between 10-30 mL/cmH 2 O [43,47]. A different technique for measuring compliance is to apply small oscillations in volume (at frequencies much higher than resting breathing frequencies) to spontaneously breathing birds, and fitting the response to a series R-I-C (resistance-inertancecompliance) model. Using this technique the 'dynamic compliance' of ducks is estimated at 7.7 mL/cmH 2 O [43]. However, this technique is very dependent on the chosen model being fitted, and relies on the overall resistance being constant throughout the breathing cycle. The differences in compliance measured with different techniques suggests that the overall compliance of the system is sensitive to body position and muscular tone.
The elastance of the airsacs is found to be approximately 1/20th of the overall elastance [49], this means that we would expect the compliance of the airsacs to be approximately 20 times higher than the compliance of the overall system (recall that compliance = 1/elastance). Because of this uncertainty in the compliances, Urushikubo et al. [30] use compliances ranging from 0.009 to 900 mL/cmH 2 O for their lumped parameters model. In this work, we find that varying the total compliance does not affect many of the qualitative dynamics of the flow. However, with low compliances (less than 100 mL/cmH 2 O), the flow through the parabronchi is sharply changing between inspiration and expiration, as shown in Fig 14. We find that compliances over 100 mL/cmH 2 O produce airflow dynamics that vary smoothly and match experimental observations [10] much better.
Ratio of airsac compliances. To determine the ratio of airsac compliances, γ, we looked at the phase of airflow through the parabronchi [10,11] and used data on the ventilation ratios of the cranial and caudal sets of airsacs [13,38,47]. Experiments found that for spontaneously breathing ducks the effective ventilation of the caudal airsacs is 56.9% of the total ventilation of all airsacs [38]. In our model, we can achieve this ventilation ratio for different γ values Fig 14. The shape of the flow q P smoothes out as C tot increases. This figure plots the flow through the parabronchi q P against time for a range of C tot , with all traces aligned such that the beginning of inspiration is at t = 0. The parameter γ = 1, other parameters are given in Table 1. Note: the transition from inspiration to expiration happens at close to the same time for C tot ! 90 mL/cmH 2 O and is shown as a black dashed line. The time of the transition for C tot = 9 mL/cmH 2 O is shown with a dark blue dot-dashed line. depending on C tot , e.g. if we choose C tot = 450 mL/cmH 2 O, we would need γ = 1.35 (the dashed line in Fig 15). However, for all compliances we found that γ > 1 (i.e., C 1 > C 2 ) for spontaneous breathing in ducks. Additionally, in relaxed (anaesthetized) conditions the ventilation ratio is 74.2% to the caudal airsacs, which is thought to be due to large increases in the compliance of the caudal airsacs [38]. To achieve this with C tot = 450 mL/cmH 2 O we would require γ = 3.81 (the dot-dash line in Fig 15), where C 1 > C 2 as would be expected.
Forcing. In all our simulations we used the sinusoidal function given in Eq (11) to model the pressure P ext outside the airsacs. The amplitude of the change in pressure due to breathing, P amp = 0.5 cmH 2 O, is chosen to match the tidal volumes [38,50] and pressures [15] seen experimentally. The period T = 3s is chosen to match the period of respiration of ducks at rest [9,13,38,50]. Using the same P ext (t) for both sets of airsacs is justified, as we can consider the thoracic-abdominal cavity (coelom) which contains all the airsacs as a single compartment with uniform pressures [51]. Experiments show that the baseline pressure in the thoracic-abdominal cavity (coelom), P c , can vary depending on physiological conditions, i.e. when crowing roosters exhibit a transient increase in P c up to 40 cmH 2 O above atmospheric pressure [52]. In this paper we choose P c = P atm for spontaneous breathing. This choice will not affect the airflow dynamics, but will alter the calculation of volumes.
Airsac volumes. The total volume of the cranial airsacs is found to be less than that of the caudal airsacs in a lot of avian species [2,53]. In ducks, the peak (end-inspiration) volumes of the caudal and cranial airsacs during spontaneous breathing were found to be 235.4mL and 221.9mL, respectively [38]. We choose V 1,res = 105.6mL and V 2,res = 103.6mL for a single side as this gives the experimental peak volumes. In this paper, with P c = P atm , the steady state volumes in the absence of changing pressures due to breathing dP ext dt ¼ 0 À Á are V 1 = V 1,res and V 2 = V 2,res . This can be calculated from Eqs (9) and (10) with P 1 = P 2 = P atm and P amp = 0. The ratio of caudal to cranial airsac ventilation increases as the ratio of airsac compliances γ = C 1 /C 2 increases, i.e. the volume of the caudal airsacs changes more than that of the cranial airsacs. The experimentally measured ratio [38] in spontaneous breathing (dashed line) and artificial ventilation (dotdashed line) ducks are shown. For different C tot values, we would require slightly different γ values to match the experimental findings. The required γ values in each case for C tot = 450 mL/cmH 2 O are shown with vertical dotted lines. All other parameters can be found in Table 1.