Computational Modelling of Cough Function and Airway Penetrant Behavior in Patients with Disorders of Laryngeal Function.

OBJECTIVE/HYPOTHESIS: Patients with laryngeal disorders often exhibit changes to cough function contributing to aspiration episodes. Two primary cough variables (peak cough flow: PCF and compression phase duration: CPD) were examined within a biomechanical model to determine their impact on characteristics that impact airway compromise. STUDY
DESIGN: Computational study.
METHODS: A Computational Fluid Dynamics (CFD) technique was used to simulate fluid flow within an upper airway model reconstructed from patient CT images. The model utilized a finite-volume numerical scheme to simulate cough-induced airflow, allowing for turbulent particle interaction, collision, and break-up. Liquid penetrants at 8 anatomical release locations were tracked during the simulated cough. Cough flow velocity was computed for a base case and four simulated cases. Airway clearance was evaluated through assessment of the fate of particles in the airway following simulated cough.
RESULTS: Peak-expiratory phase resulted in very high airway velocities for all simulated cases modelled. The highest velocity predicted was 49.96 m/s, 88 m/s, and 117 m/s for Cases 1 and 3, Base case, and Cases 2 and 4 respectively. In the base case, 25% of the penetrants cleared the laryngeal airway. The highest percentage (50%) of penetrants clearing the laryngeal airway are observed in Case 2 (with -40% CPD, +40% PCF), while only 12.5% cleared in Case 3 (with +40% CPD, -40% PCF). The proportion that cleared in Cases 1 and 4 was 37.5%.
CONCLUSION: Airway modelling may be beneficial to the study of aspiration in patients with impaired cough function including those with upper airway and neurological diseases. It can be used to enhance understanding of cough flow dynamics within the airway and to inform strategies for treatment with "cough-assist devices" or devices to improve cough strength. LEVEL OF EVIDENCE: N/A.

INTRODUCTION

Cough is an airway protective function serving both preventative and corrective roles in pulmonary health. Physiologically, cough is a coordinated series of respiratory, laryngeal, and pharyngeal muscle activity. It can be identified and characterized perceptually, acoustically, and through direct sampling of airway pressures and flows.1, 2, 3, 4, 5 The cough aerodynamic sequence of inspiration, compression, and expulsion protects the lungs through removal of secretions and foreign material from the airway. Patients with laryngeal/airway disorders often exhibit weak and slow cough, which may result in aspiration. A weak and slow cough response is impacted by two parameters: the compression phase duration (CPD), which represents the period of full glottal closure immediately prior to cough and the peak‐cough flow (PCF), an indirect measure of cough strength reflecting the velocity of expiratory airflow.6, 7

Common treatment targets for patients with reduced airway clearance focus on improving CPD and PCF. Some of these include: postural changes, incentive spirometry, expiratory muscle strength training, percussion and vibration, manually assisted cough, huffing, and active breathing techniques. Many of these strategies include maneuvers resulting in forced expiration, creating greater “shearing” forces within the airway for elimination of foreign material.8

Recent advances in computational techniques and medical imaging have enabled the integration of engineering principles with medicine in the fast‐developing field of biomechanical modelling. Biomechanical models derive from the laws of physics and provide cost effective, non‐invasive techniques for simulation of biological processes and systems in a manner that is often difficult to accomplish in a clinical setting.

The purpose of this study was to develop a preliminary biomechanical model to simulate cough strength and aspiration. To build the model, calculation of the spatio‐temporal flow‐particulate behavior was developed through application of transient flow rates by manipulation of CPD and PCF. It was hypothesized that the model would be able to track the trajectory of airway penetrants, and be used to quantify aspiration secondary to changes in CPD and PCF.

MATERIALS AND METHODS

Biomechanical Model

The representative 3D airway geometry was constructed from computed tomography (CT) scan images. The image is segmented in order to extract the 3D geometry of the upper airway associated with a cough event, by use of MIMICS (Materialise; Leuven, Belgium), a widely applied medical image segmentation software. After segmentation, the reconstructed geometry is divided into small finite volumes or cells called meshes by use of 3‐MATIC (Materialise; Leuven, Belgium) commercial software code. Due to the irregular airway geometry, an unstructured mesh system was used to allow mapping of the exact shape. A combination of triangular (surface mesh) and tetrahedral (volume mesh) elements were used to capture geometrical features. The number of faces, vertices, and mesh cells used in the study were 837,223, 2,425,705, and 810,665, respectively. A grid‐independence test was first carried out using three different grid sizes (∼ 600,000, 800,000, and 1,000,000 cells) on order to ensure accuracy of the results. The computed results were found to be essentially independent of the grid size beyond 800,000 cells, and this grid structure was used for all results presented in the paper.

The equations governing fluid flow and particulate dynamics were integrated over the finite volumes generated by means of STARCCM+ commercial software (CD‐adapco; Melville, New York, U.S.A.). The velocity was first calculated by solving two conservation equations derived from physical laws, namely, conservation of mass (mass of the fluid remains unchanged as it flows through the airway) and conservation of momentum or Newton's second law of motion for fluid element.

Mass conservation:

Momentum equation: where, , represent time‐averaged air flow velocities in x, y, and z directions, is dynamic viscosity, and is turbulent (eddy) viscosity. The term represents the momentum exchange between the dispersed phase (particulate) and carrier gas (air) and includes parts of the stress tensor not appearing directly in the diffusion and pressure gradient terms of the fluid flow equation. The dynamics of penetrant particle in the cough flow are governed by the sum of all forces acting on the particle using Newton's second law of motion. where is the particle velocity, is mass of the penetrant, is Stokes drag force acting on the surface of the penetrant, is pressure gradient force developed by undisturbed flow, is added‐mass force created by displaced fluid, and is Buoyancy force. The left‐hand side of the equation represents the inertial force acting on the penetrant due to its acceleration.

Gas density and temperature are assumed to be constant. The laminar‐turbulent transitional nature of cough flow is accounted for by use of appropriate turbulence model for calculation of the surface force in the momentum equation.9, 10, 11 Three major approaches are typically used for simulating such transitional flows: direct numerical simulation (DNS), large‐eddy simulation (LES), and Reynolds‐averaged Navier–Stokes (RANS) models.12 The k‐ω shear stress transport turbulence (SST) model in which k is the turbulent kinetic energy, and ω is the specific dissipation rate has been adopted in this study because DNS or LES turbulence models are computationally intensive. In addition, the k‐ω SST model has been successfully used for transitional flows, and represents an adequate compromise between computational intensity and accuracy.13

The locations of individual particles are tracked while being transported through the fluid in the airway, a method known as Lagrangian approach.14 The particle dynamic equation is augmented by a sub‐model allowing for possible droplet break‐up, collision and coalescence. The underlying mechanism is that the surface tension and viscous forces inside the droplet resist its deformation.15, 16

The governing equations are solved subject to appropriate boundary and initial conditions relevant to cough. Figure 1 shows the airway geometry indicating the inlet and outlet boundaries and the other upper airway landmarks used. The outlet is defined at the top surface of the airway geometry where zero gage pressure is imposed relative to atmospheric pressure. This approach expresses the fact that the mouth is typically open during cough. A voluntary cough waveform (see Fig. 2) obtained from a human subject was applied at the inlet defined at the two (bifurcation) interfaces between the trachea and the lungs.

Figure 1

Airway geometry indicating inlet and outlet boundary surfaces and other landmarks

Figure 2

Cough waveform obtained from human subject and used as input data for the simulation

Parametric Study

A key objective of this study was to assess the effect of two cough parameters: compression phase duration (CPD) and peak cough flow (PCF), which determine the quality of cough strength. Although a typical voluntary cough has three phases: inspiratory, compression, and expiratory, most patients do not inspire prior to cough when a laryngeal penetrant is present.17 Therefore, for the purpose of this study, the simulation was started from the compression phase of cough as shown in Figure 2.

Liquid droplet penetrants were considered to have the properties of a “thin liquid,” e.g., water. The droplet penetrants were introduced at four different landmarks (levels) along the laryngeal section of the airway as shown in Figure 1. The laryngeal section of the airway possesses abundant sensory innervation, which can activate to produce cough.18 These landmarks represented the levels of the laryngeal vestibule, above the laryngeal vestibule, vocal folds and above the vocal folds where cough receptors are located.19 For the purpose of this simulation study, a single droplet penetrant with 4 mm diameter was introduced at each level of both the anterior and posterior walls of the airway. The stickiness of the airway wall due to the presence of mucus was represented in the biomechanical model through a wall restitution parameter, which typically ranges from 0 (complete stickiness) to 1.0 (no stickiness). For this study we adopted a value of 0.1 in order to allow for a possible mucus‐induced high degree of stickiness.20

In order to establish the effects on the penetrants of the key parameters considered (airflow velocity distribution, penetrant trajectory, and penetrant break‐up pattern), a simulated data set was first generated. A number of studies have shown the effect of PCF and CPD in different population groups. Specifically, Pitts et al. (2008) have shown that the penetration/aspiration groups with Parkinson's disease produce significantly longer CPD time.21 Also, Kulnik et al. (2016) has found that higher PCF lowers the risk of pneumonia in acute stroke.22 The percentage change associated with PCF and CPD in these studies vary significantly according to the population considered. Therefore, our study is designed to represent all possibilities associated with the key parameters considered. The base case represents an impaired voluntary cough from a human subject. Four cases (Cases 1 to 4), were then generated by changing both the CPD and PCF ± 40%.23 Case 1 demonstrates both CPD and PCF as decreased by 40% compared to the base case, Case 2 demonstrates the CPD reduced by 40% with an increase of PCF of 40%. Case 3 and Case 4 simulate a 40% increase to CPD with a corresponding PCF decrease of 40% in Case 3 and increase of 40% in Case 4.

RESULTS

The results of this study can be evaluated across three specific parameters; the velocity distribution, the penetrant trajectory and the pattern of break‐up of the penetrant once exposed to the cough event.

Velocity Comparison

Predicted transient variation of velocity at the level of the vocal folds for the base case and four simulated cases (Case 1‐ 4) are presented in Figure 3. In all cases, the highest velocity was observed at the peak expiratory flow. The highest velocity predicted for these cases was 49.96 m/s, 88 m/s, and 117 m/s (Cases 1 and 3, Base case, and Cases 2 and 4) respectively. As demonstrated in Figure 3, the predicted transient velocity pattern followed the same profile as the flow rate.

Figure 3

Predicted velocity magnitude at the level of the vocal folds

Penetrant Trajectory

Figure 4 shows the trajectories of penetrants released at location P1 for the base case and four simulated cases considered in the study. P1 is located on the posterior side of the airway above the laryngeal vestibule. The real‐time duration of the simulation (Ts) in each case from the instant of release is also shown in Figure 4. The final location of the penetrant is indicated as FP in each figure.

Figure 4

Predicted trajectories of liquid penetrants released at P1 above the larynx

Results indicate that the penetrant was retained in the airway in all cases except Case 2. For Cases 3 and 4, the prolonged CPD resulted in the penetrant falling further down the airway before the expiratory phase. If the penetrant travelled to the lung‐end of the airway it was considered “aspirated.”

The trajectory of penetrants considered in all cases is summarized in Table 1. The three major outcomes predicted in the study are identified by the letters E (for particles that have escaped from the airway to the oral cavity), R for particles retained in the airway, and A for those aspirated to the lungs. We define three indices that may quantify aspiration, retention, and clearance of penetrants. For our model, we have assumed that cough is considered a success if penetrants are cleared from the laryngeal airway. Likewise cough is considered as “in need of a sequential cough” if penetrants are retained. The last parameter defines a cough “failure” if penetrants are aspirated.

Table 1

Penetrant trajectory outcomes and evaluated indices (Key: R‐Retained in airway, E‐Escaped to oral cavity, and A‐Aspirated to lungs)

Penetrant	Base case	Case 1	Case 2	Case 3	Case 4
P1	R	R	E	R	A
P2	R	E	R	R	E
P3	E	R	E	R	E
P4	R	A	R	A	A
P5	E	E	E	E	E
P6	R	A	R	A	A
P7	A	E	E	A	A
P8	A	A	R	A	A
success index = number of escaped penetrantstotal number of penetrants	0.25	0.375	0.5	0.125	0.375
need of sequential cough (NSC) index = number of retained penetrantstotal number of penetrants	0.5	0.25	0.5	0.375	0
failure index = 1−(success+NSC)	0.25	0.375	0	0.5	0.625

In the base case, 25% of the penetrants escaped into the laryngeal airway. The highest percentage (50%) of penetrants escaping into the airway were observed in Case 2 (with −40% CPD, +40% PCF), while only 12.5% escaped in Case 3 (with +40% CPD, −40% PCF). The proportion that escaped in Cases 1 and 4 was 37.5%.

The results of the penetrant trajectories provided in Table 1 reveal that penetrant P5 was successfully cleared from the airway in all simulated cases, while P4, P6, and P8 (introduced from the anterior airway) were not. Figure 5 shows the difference in penetrant trajectory for the same level when introduced anteriorly versus posteriorly.

Figure 5

Penetrant trajectories for two penetrants (P6: anterior and P5: posterior)

Penetrant Break‐up Pattern

An example of a penetrant break‐up event within the airway is presented for location P5 in Case 3 (Fig. 6). The penetrant released at location P5 was predicted to creep down the wall of the airway due to gravity during the compression phase. Once the expiratory phase is initiated, the penetrant is disturbed by the exchange of momentum with the increasing airflow, resulting in disintegration of the droplet into small (child) droplets. The number of child droplets created during this “break‐up” event depends on several factors: 1) the droplet resident location, 2) the relative speed of the air flow, 3) the initial droplet size, and 4) the surface tension. The results presented in Figure 6 were obtained using the best estimates of these parameters and are considered to be only illustrative of droplet behaviour at this stage of the study.

Figure 6

Predicted fate of liquid penetrant released from location P5 for Case showing cough‐induced disintegration into small droplets

Figure 7 models 6 instances where the droplet break‐up event is visualized using representative spherical particles for the penetrant P5 considered in Figure 6, based on the simulation. The real‐time duration after cough initiation is shown in time steps of 0.001 s from 1.180 s to 1.18 5s. The penetrant is initially represented as one droplet, which upon subsequent interaction with the flow velocities of about 40 m/s, breaks up into small droplets. These child droplets then travel up the airway towards the mouth‐end. This process is illustrated in Figures 7(b)‐7(f). It should be noted that the child droplet trajectories and velocities are different from each other as they travel along the airway.

Figure 7

A droplet break‐up event visualized at 6 instances

The mean diameter of the penetrant changes during its residence within the airway for the time steps corresponding to those in Figure 7. At the first breakup event, the creation of child‐droplets reduces the mean diameter of the liquid penetrant from initial 4 mm to 0.91 mm. The mean diameter is subsequently reduced to 0.38 mm as a result of the second break‐up event at 1.183 s. While travelling through the airway, some of these child‐droplets may also collide and coalesce with each other due to change in the flow characteristics to form new droplets. Such coalescence may increase the mean diameter of the child droplet as predicted at 1.185 s, and the new mean diameter is now 0.62 mm.

DISCUSSION

This study investigated the effect of CPD and PCF variations in voluntary cough function using biomechanical modelling within the context of a 3D airway geometry reconstructed from CT‐scan images of human airways. The results of this study elucidate cough function and penetrant behavior across three informative parameters: 1) airflow velocity distribution, 2) penetrant trajectory, and 3) break‐up pattern. The findings suggest that specific modifications to the expected cough dynamics may promote airway clearance.

Airflow Velocity

Specifically, the analysis of airflow velocity revealed that the modification in PCF resulted in substantial change in airway velocities for all cases modelled. This is consistent with previous studies and confirms the value of peak cough flow measures as a surrogate for velocity within an airway clearance event.24, 25 The largest velocity magnitude is recorded in Case 2 and Case 4 whereas the smallest is recorded in Case 1 and Case 3. The transient velocity variation in Case 1 and Case 2 is advanced in time compared to other cases due to short CPD. Therefore, the Case 2 demonstrates a strong fast cough response whereas Case 3 can be considered as a weak, slow cough. The outcome of velocity profiles show that CPD primarily manipulates the pace of cough (fast/slow) while PCF contributes to the power of the cough (strong/weak).

The trajectory study showed that Case 2 cleared more penetrants than any other cases. Recall that in Case 2, the CPD is decreased by 40% from the base case, and PCF is increased by 40%. The peculiarity of Case 2 appears due to the fact that the combined effects of these changes to CPD and PCF cause rapid flow acceleration after the short CPD duration, resulting in favorable conditions for the penetrant to break up before escaping from the mouth‐end of the airway. In contrast, when the CPD was increased by 40% in Cases 3 and 4, the prolonged CPD resulted in the penetrant falling further down the airway before the expiratory phase. In addition, our data revealed that penetrants originating from the anterior side of the airway are likely to fall deeper into the airway, presenting a greater challenge for airway clearance. A comparison of an anterior and a posterior penetrant at the same level in the airway showed the fact that the airway geometrical shape does not obstruct the path of the penetrants originating from anterior side. These data are novel and may provide direction to clinicians tracking events using dynamic video measurement. Consequently, the parameters defined by the simulated cases considered in this study may be utilized to direct therapeutic cough‐assist applications. The findings of the study also suggest that cough‐assist applications may need to effectively increase PCF and reduce CPD in order to maximize cough effectiveness as a defensive mechanism against aspiration.

Penetrant Break‐Up Pattern

Finally, our study also comprehensively modelled penetrant break‐up. During a break‐up event, the penetrant diameter significantly drops when small child‐droplets are created. In addition, the collision of liquid droplets together may enhance coalescence and increase the mean diameter of the droplets residing within the airway. By comparing the mean diameter, we may be able to predict the number of child droplets created in a specific break‐up event. Previous studies measured aerosol droplet sizes leaving the mouth, but not within the airway.26 Our findings shed light on a previously unrecognized issue within cough evaluation and management interventions. The role of droplet collision, coalescence and re‐formation as a risk marker, and the impact of scatter for micro‐aspiration events may need further consideration.

The data presented from our modelling efforts are both novel and informative to cough‐assist evaluation and interventions. Modelling can assist in specifying different classes of features previously not recognized in contributing to event outcomes. Models like ours presented here can help us to better understand and possibly predict the behaviour of cough function in conditions that cannot easily be experimentally reproduced. Despite this, the diversity of system conditions, the availability of system‐specific information and availability of patient specific data may limit outputs from any model and imply the need for validation and further adaptive models. One stiff assumption in our study is the use of static CT images to reconstruct the 3D airway geometry which thus has not considered the structural changes of the laryngeal wall. The use of optical flow‐based motion estimation, which is derived from 4D CT image data set will enable representation of the dynamic nature of the airway wall in subsequent studies. Like other modelling techniques, our efforts here must be further tailored not only to make the best use of the information available but also to answer additional questions such as; how does cough bolus or material type interact with penetrant trajectory? How will anatomical deviations (i.e., cervical spondylosis) impact the velocity and tracking of the penetrants? Likewise, some of the parameters altered in our models were based upon clinical knowledge of cough and patient cough experience. As a result some bias and inaccuracy may be present in those decisions. Similarly we have not considered the interactive events of cough and swallowing in airway clearance and penetrant trajectory. Although many questions remain, this study presents “proof of concept” in the development of a cough efficiency computational model.

CONCLUSION

Airway modelling may be beneficial to the study of aspiration in patients with impaired cough function including those with upper airway and neurological diseases. It can be used to enhance the understanding of cough flow dynamics within the airway and to inform strategies for treatment with “cough‐assist devices” or devices to improve cough strength. The proposed biomechanical model is the first of its kind to study cough‐penetrant dynamics within the human airway. As such it presents a “proof‐of‐concept” in the biomechanical modelling of human cough.