Modeling Lung Derecruitment in VILI Due to Fluid-Occlusion: The Role of Emergent Behavior

Ventilator-induced lung injury (VILI) is driven by the processes of volutrauma and atelectrauma, which can act synergistically to compromise the blood-gas barrier. We have postulated that this synergy arises through a rich-get-richer mechanism whereby atelectrauma causes holes to form in the blood-gas barrier while concomitant volutrauma causes susceptible holes to progressively enlarge as VILI worsens. We previously developed an analytical model based on this idea that accurately predicts the progressive increases in lung elastance seen immediately following a recruitment maneuver as VILI progresses over the course of hours. In the present study we extend this model to account for the rate of change of elastance, due to closure of lung units, in the minutes following a recruitment maneuver. We found that the distribution of unit closing velocities throughout the lung can be described by a power law with an exponent of −2 that matches previously published power laws associated with the dynamics of lung recruitment. Our model thus reveals lung collapse as an example of emergent complex behavior and links the dynamics of altered function in the injured lung to structural damage in a way that explains the mechanisms of injury progression arising from the ongoing stresses and strains applied by mechanical ventilation.

Introduction

Mechanical ventilation plays a major role in critical care, providing life support to patients in respiratory failure. However, because of the stresses and strains it imposes on the tissues of the lung, mechanical ventilation can also cause ventilator-induced lung injury (VILI) (Dos Santos and Slutsky, 2000; Gattinoni et al., 2003; Slutsky and Ranieri, 2013). The two principle injury mechanisms of VILI are over-distension of the tissues that gives rise to volutrauma, and cyclic recruitment and derecruitment of lung units that gives rise to atelectrauma (Parker et al., 1993; Dreyfuss and Saumon, 1998). Volutrauma and atelectrauma are both involved in compromising the blood-bas barrier of the lung, which allows plasma-derived fluid and proteins to leak into the airspaces. This disrupts surfactant function, which then decreases lung compliance and thus increases the tissue stresses associated with mechanical ventilation, making the leak worse. At the same time, the lung becomes progressively more derecruited, which places the remaining open lung regions at greater risk of over-distension (Gattinoni and Pesenti, 2006; Gattinoni et al., 2016). Patients with pre-existing lung injury, particularly acute respiratory system distress syndrome (ARDS), are highly susceptible to VILI (Slutsky and Ranieri, 2000; Dreyfuss and Hubmayr, 2016). Various ventilation strategies have been proposed to manage ARDS and prevent VILI, but once underway it is difficult to break out of the vicious cycle of progressing injury (Slutsky and Ranieri, 2013).

We have shown in a number of studies in mice (Seah et al., 2011; Smith et al., 2013, 2017) that the nature and degree of lung injury is reflected in what we term the derecruitability of the lung. Derecruitability is measured by applying a recruitment maneuver (deep inflation) to maximally open closed airspaces followed immediately by a period (typically 3 min) of mechanical ventilation at a designated tidal volume (V_t) and level of positive end-expiratory pressure (PEEP). Lung stiffness (H) rises progressively during the ventilation period to a degree that is modest (perhaps 20%) in a normal lung, but is dramatically enhanced in the injured lung to an extent that reflects the degree of injury (Smith et al., 2013). Some of the post-recruitment increase in H in the normal lung likely reflects changes in surface tension arising from the dynamics of surfactant at the air-liquid interface and viscoelastic adaptation of tissue stress, but the highly exaggerated increases in H seen in the injured lung are largely the result of progressive derecruitment of airspaces (Massa et al., 2008; Smith et al., 2017; Hamlington et al., 2018b).

We denote as D_rate the mean rate of increase in H over the 3 min following recruitment, while H₁ is the first measurement of H made immediately after recruitment (Smith and Bates, 2013; Smith et al., 2013). Both D_rate and H₁ increase with the severity of lung injury, but they do not mirror each other exactly. D_rate quantifies derecruitment dynamics that manifest over a timescale of minutes, probably reflecting instabilities in the fluid layer lining the small airways and alveoli (Smith et al., 2013). H₁, on the other hand, quantifies either much more rapid derecruitment phenomena or an inability to recruit at all, more likely explicable on the basis of rapid alveolar flooding by accumulated edema fluid and/or unstable collapse of lung units caused by high surface tension in the air-liquid interface (Smith et al., 2013). D_rate and H₁ together thus comprise a pair of sensitive biomarkers of lung injury (Smith et al., 2013).

We have previously shown in mice that, as VILI develops, the time-course of H₁ can be accurately accounted for in terms of the accumulated leak of plasma-derived material through holes in the blood-gas barrier (Mori et al., 2018). These holes appear to be generated by both atelectrauma and volutrauma acting synergistically in a rich-get-richer fashion (Hamlington et al., 2018a). The rate of accumulation of edema in the airspaces, and thus the displacement of air that results, is a function of the number and sizes of these holes, both of which increase with time during injurious mechanical ventilation. In the present study, we extend this analysis to D_rate in order to establish a consistent theory of derecruitability of the injured lung that links both H₁ and D_rate to the biophysical processes taking place during the development of VILI.

Materials and Methods

Model Development

In the following model development we take a similar analytical approach to our previous study of the dependence of H₁ on degree of lung injury, this time focusing on the complementary parameter D_rate. To account for the progressive increase in D_rate that accompanies the development of VILI, it is necessary to model events taking place over two different scales of time. The longer time scale, on the order to minutes to hours, concerns the progression of VILI itself. The shorter time scale, on the order of seconds to minutes, involves the dynamic recruitment and derecruitment of alveoli and small airways that are reflected in our periodic assessments of VILI. In order to avoid ambiguity between these two time scales in the following mathematical development, events taking place over the shorter time scale (seconds to a few minutes) will be referenced to time denoted by τ while events taking place over the longer time scale (many minutes to hours) will be referenced to time denoted by t.

We first consider the short time-scale assessment of VILI, which is based on the computational model of recruitment/derecruitment dynamics proposed by Bates and Irvin (2002) consisting of a parallel distribution of lung units. Each unit has its own airway to which is applied the same pressure. The units open and close according to their randomly assigned critical opening and closing airway pressures (P_o and P_c, respectively). These pressures are drawn from probability distribution f_Po(P_o) and f_Pc(P_c), respectively. However, a unit does not necessarily open or close as soon as it traverses its particular value of either P_o or P_c. Rather, a latency is built into the process by associating with each unit a virtual trajectory variable x that can assume any value between 0 (corresponding to the unit being definitively closed) and 1 (definitively open). An open unit remains open until x=0, while a closed unit remains closed until x=1. Movement of x along the virtual trajectory occurs to the right (i.e., toward increasing values of x) when P > P_o, and to the left when P < P_c, with the rate of change x being proportional to the amount by which P either exceeds P_o or is less than P_c. The constants of proportionality that determine the velocity of x rightward and leftward are s₀ and s_c, respectively. That is,

$$\frac{dx}{d\mathrm{\tau }}\mathrm{\_}=\left(P-P_o\right)s_o;\text{for}P\ge P_o$$

$$=\left(P-P_c\right)s_c;\text{for}P\le P_c$$

$$=0;\text{for}{P}_c<P<P_o(1)$$

where 0≤x <  1. The velocity constants s₀ and s_c for each unit are drawn from probability distributions f_so(s₀) and f_sc(s_c), respectively.

Following the original model (Bates and Irvin, 2002), we assume that P_o and P_c are distributed according to Gaussians with variances σ₀ and σ_c, and means μ₀ and μ_c, respectively. In accord with the findings of Massa et al. (2008) in acid-injured mice, we set the two variances to be equal (i.e., σ₀ = σ_c = σ) and the two means separated by δP such that the opening pressures are greater than the closing pressures (i.e., μ₀ = μ_c + δP)

$$f_{P_c}\left(P_c\right)=\frac{1}{\mathrm{\sigma }\sqrt{2\mathrm{\pi }}}{e}^{-\frac{1}{2}{\left(\frac{P_c-{\mathrm{\mu }}_c}{\mathrm{\sigma }}\right)}^2}(2)$$

and

$$f_{P_o}\left(P_o\right)=\frac{1}{\mathrm{\sigma }\sqrt{2\mathrm{\pi }}}{e}^{-\frac{1}{2}{\left(\frac{P_o-{\mathrm{\mu }}_c-\mathrm{\delta }P}{\mathrm{\sigma }}\right)}^2}(3)$$

We assume that σ and δP do not change with increasing injury, so the progression of VILI is reflected entirely in the way that μ_c increases with time.

The velocity distributions f_so(s₀) and f_sc(s_c) were originally specified on purely empirical grounds as conforming to hyperbolic functions (Bates and Irvin, 2002), motivated by the ubiquitous appearance of such functions in naturally occurring complex systems (West, 1990; Robertson et al., 2000; West et al., 2000). Solving the model analytically using such functions is not possible, however, because hyperbolas are not finite-integrable over s ∈ (0,∞), but finite integrals can be obtained by having the exponent of the power law be different from unity and by specifying a small positive lower limit, s_m, for integration. Accordingly, we set

$$f_{s_c}(s_c)=\frac{\mathrm{\alpha }-1}{s_m}{\left(\frac{s_c}{s_m}\right)}^{-\mathrm{\alpha }}(4)$$

where 0 < s_m < s_c < ∞ and α > 1.

When modeling how the fraction of closed lung, F_closed, changes during a derecruitability test, we avoid having to consider f_so(s₀) specifically by assuming, as did Albert et al. (2009), that it is only derecruitment that takes place to any significant extent during the 3 min following a recruitment maneuver. That is, we assume that the fully recruited lung derecruits progressively during these 3 min as if exposed to an effective pressure, P_derecruit, that is fixed throughout the breath and that reflects the mean airway pressure applied by the mechanical ventilator. P_derecruit is determined by the applied level of PEEP plus some positive ΔP that depends on the ventilation regimen that is applied on top of the PEEP. As the lung becomes progressively more injured, and thus stiffer, mean airway pressure increases. This pressure increase is offset to some extent, however, by the compliance of the gas in the piston and connecting tubing of the flexiVent ventilator (Scireq, Montreal, Canada) that was used to collect the experimental data. While the effect of this compliance on calculation of H was compensated for digitally (Schuessler and Bates, 1995), it still absorbed an increasing fraction of the pressure increases as the lung stiffened (Smith et al., 2013). Mean airway pressure also increases somewhat over the course of a single derecruitability test. However, in order to obtain an analytical solution to our model we assumed that P_derecruit remains unchanged over each of the 3 min derecruitment maneuvers. This is clearly not precisely the case, but we take it to be a reasonable first-order approximation since tidal volume is quite low during the tests.

For a unit to have the possibility of closing during a derecruitability test it must therefore satisfy the condition P_c > P_derecruit. Whether or not a particular unit actually closes during the test, however, depends on whether it moves along its virtual trajectory rapidly enough for x to reach a value of 0 before the 3 min test is complete and that there is no significant retrograde motion of x during the inspiratory portions of the test. If the duration of the test were to continue indefinitely, the asymptotic value to which F_closed would tend is that which satisfies the above condition, which is

$$F_{\text{closed}}\left(\mathrm{\infty }\right)=\underset{P_{\text{derecruit}}}{\overset{\mathrm{\infty }}{\int }}\left[\frac{1}{\mathrm{\sigma }\sqrt{2\mathrm{\pi }}}{e}^{-\frac{1}{2}{\left(\frac{P_c-{\mathrm{\mu }}_c}{\mathrm{\sigma }}\right)}^2}\right]d{P}_c(5)$$

This situation is never actually achieved, of course, because the derecruitability test is terminated after 3 min.

We assume that the recruitment maneuver establishes an initial condition for which all units open and are positioned at x=1 along their respective virtual trajectories. As derecruitment proceeds during the test, the various units in the model close sequentially. The time, τ, at which a given unit closes is determined by how quickly it traverses the length of its virtual trajectory from x=1 to x=0, which is given by

$$\mathrm{\tau }=\frac{1}{s_c\left|P_{\text{derecruit}}-P_c\right|}(6)$$

which means that the units that have closed by time τ are those that have closing pressures satisfying $\frac{1}{\mathrm{\tau }\left(P_c-P\right)}<s_c$. Assuming, as did Albert et al. (2009), that f_Pc(P_c) and f_sc(s_c) are statistically independent, F_closed(P,τ) is given by

$$\begin{array}{cc}{F}_{\text{closed}}(P,\mathrm{\tau })=\underset{P}{\overset{\mathrm{\infty }}{\int }}\underset{{\left[\left(P_c-P\right)\mathrm{\tau }\right]}^{-1}}{\overset{\mathrm{\infty }}{\int }}\frac{1}{\mathrm{\sigma }\sqrt{2\mathrm{\pi }}}{e}^{-\frac{1}{2}{\left(\frac{P_c-{\mathrm{\mu }}_c}{\mathrm{\sigma }}\right)}^2}\hfill & \\ \cdot \frac{\mathrm{\alpha }-1}{s_m}{\left(\frac{s_c}{s_m}\right)}^{-\mathrm{\alpha }}d{P}_{c}d{s}_c\hfill & \end{array}(7)$$

A detailed step-by-step solution of Eq.7 is presented in the Supplementary Material.

Solving Eq. 7 we find that α = 2 and Eq. 4 of thus becomes

$$f_{s_c}(s_c)=\frac{1}{s_m}{\left(\frac{s_c}{s_m}\right)}^{-2}(8)$$

Massa et al. (2008) found σ = 3 in mice with acute lung injury caused by hydrochloric acid instillation. Smith et al. (2013) found μ_c(t) = kt, where k is a constant, in injuriously ventilated mice. Since PEEP is zero and H is measured using small-amplitude perturbations in lung volume, we assume that P=0 to a first approximation (Mori et al., 2018). D_rate can thus be written as:

$$D_{\text{rate}}\left(t\right)=\frac{H_1\left(t\right)s_m}{\sqrt{2\mathrm{\pi }}}\left[\sqrt{\frac{\mathrm{\pi }}{2}}3kt\left(erf\left(\frac{kt}{3\sqrt{2}}\right)+1\right)+9e^{-\frac{1}{2}{\left(\frac{kt}{3}\right)}^2}\right](9)$$

From our previous proposed model (Mori et al., 2018) we can relate Eq. 9 to H₁ when t is sufficiently large:

$$\underset{t\to \mathrm{\infty }}{lim}{D}_{\text{rate}}\left(t\right)=3k{s}_m\frac{1}{\sqrt{a}}\left(H_1\left(t\right)-\frac{H_0}{2}\right)(10)$$

Experimental Data

We tested the predictions of the above model against a set of data collected in mice and published previously (Smith et al., 2013). A detailed description of the experimental method used, as well as animals ethics approval, are given in this previous publication. Briefly, healthy 8- to 10-week-old female BALB/c mice (18.3–24.4 g) were anesthetized with 90 mg/kg intraperitoneal sodium pentobarbital and then connected to a computer-controlled mechanical ventilator (flexiVent, Scireq, Montreal, Canada). Following paralysis (0.5 ml/kg IP injection of pancuronium bromide) animals received a 21 min protocol starting with 16.5 min of injurious ventilation delivered at zero PEEP with a very large V_t (1.0, 1.1, 1.2, and 1.3 ml in 4 different groups of animals) followed by a derecruitability test (4.5 min protocol consisting of a deep inflation followed by ventilation at zero PEEP with V_t = 0.25 ml during which elastance was measured every 20 s). This 21 min protocol was repeated continuously for 4 h or until the animal died. The way that the lungs derecruit over time during the derecruitment test is reflected in the way that H increases with time (Allen and Bates, 2004). Finally, the resulting estimations of D_rate from the derecruitability tests were fitted to Eq. 15, and the 95% joint confidence regions for s_m and k were determined as described by Mori et al. (2018).

It should be pointed out that we also measured parameters related to the resistance of both the conducting airways and the respiratory system tissues in these mice. We focus here exclusively on the elastance parameter H in our model because it scales inversely with the fraction of lung that remains open, provided intrinsic tissue stiffness does not change, thus serving as a convenient biomarker of derecruitment. Lung tissue resistance scales similarly to elastance (Smith et al., 2020) and therefore adds no new information, so we did not include it in our model. Airway resistance tends to correlate less well with lung derecruitment because the proximal branches of the airway tree, which contribute much to airway resistance, mostly remain patent even when significant portions of the lung periphery derecruit (Smith et al., 2020).

Model Fitting

The model described by Eq. 9 was fit to the elastance data versus time collected within each derecruitment maneuver with MATLAB (The Mathworks, Natick, MA, United States) using a non-linear least squares method (trust-region algorithm) for each animal individually in all 4 groups, resulting in a pair of values of s_m and k for each animal. The maximum number of evaluations and iterations were set to 600 and 400, respectively, and the termination tolerance value of the cost function (mean squared residual) was 10^–6. Finally, the 95% joint confidence regions for s_m and k were determined as described by Mori et al. (2018).

Results

Figure 1 shows the fits of Eq. 9 to experimental data from four groups of over-ventilated mice (V_t = 1.0, 1.1, 1.2, and 1.3 ml, with zero PEEP), together with the experimental data. The model fits pass within one standard error either side of the data points in most cases, and are very close to these ranges in the other cases. The R² values for all fits were greater than 0.85. This demonstrates that the model accurately accounts for the way that injury, as reflected by H, accelerates over time and how this acceleration increases with increasing V_t.

FIGURE 1

Figure 1. D_rate (average ± SE) versus time as VILI developed in the 4 groups of mice (symbols) together with the model fits to Eq. 9 (lines). Black: V_t = 1.0 ml (RMSE = 0.35); Red: V_t = 1.1 ml (RMSE = 0.64); Green: V_t = 1.2 ml (RMSE = 0.36); Blue: V_t = 1.3 ml (RMSE = 0.65). The data were collected in a previous study (Smith et al., 2013) (RMSE – root mean squared residual.).

Figure 2 shows the two model parameters in Eq. 9, s_m and k, versus V_t. There was no dependence of s_m on V_t by ANOVA (p = 0.11, Figure 2A), whereas k was significantly greater for the lungs exposed to the two higher V_t compared to the two lower V_t (p < 0.001, Figure 2B). Normality of the distributions of all parameter groups was supported by the Shapiro-Wilk test (p > 0.05). The Bonferroni post hoc test was performed to assess inter-groups differences for k. All statistical tests were performed in MATLAB 2019b (The Mathworks, Natick, MA, United States).

FIGURE 2

Figure 2. Mean ± 2 SEM for the two key model parameters s_m (A) and k (B) and for the 4 groups, 1.0 ml (n = 6), 1.1 ml (n = 6), 1.2 ml (n = 5) and 1.3 ml (n = 4). Statistically significant differences between groups are indicated by: *P < 0.05), **P < 0.01.

The joint 95% confidence intervals show clean separation between the parameter pairs for the four groups of mice (Figure 3).

FIGURE 3

Figure 3. Joint paired 95% confidence regions for the model parameters. The optimum parameter values (those that minimize the root mean squared difference between measured and modeled D_rate) are located at the centers of the ellipses. Black: V_t = 1.0 ml; Red: V_t = 1.1 ml; Green: V_t = 1.2 ml; Blue: V_t = 1.3 ml.

Figure 4 shows plots of D_rate(t) versus H₁ for all animals studied. Also shown are the model predictions of D_rate(t) from Eq. 9 versus H₁ calculated according to Mori et al. (2018). Lastly, Figure 5 shows, in a Bland-Altman plot, the difference between modeled and measured D_rate values as a function of H₁.

FIGURE 4

Figure 4. Rate of change of lung stiffness (D_rate) versus initial measurement of post-recruitment stiffness (H_1) measured in the 4 groups of mice (symbols) together with the model fits (lines). Black: V_t = 1.0 ml; Red: V_t = 1.1 ml; Green: V_t = 1.2 ml; Blue: V_t = 1.3 ml.

FIGURE 5

Figure 5. Difference between modeled and measured D_rate as a function of H_1 (open dots). Mean is represented by the full line whereas 95% interval (± 1.96 SD) is represented by dashed lines.

Discussion

We have developed an analytical model of the time-course of derecruitment in the injured lung that mimics observations of how lung derecruitability evolves during the development of VILI. This analytical model is based on our previous numerical approach of ascribing recruitment and derecruitment dynamics to the behavior of a virtual trajectory associated with various lung units (Bates and Irvin, 2002; Massa et al., 2008). This model uses a similar analytic approach to our recent study that examined how H measured immediately after a recruitment maneuver (H₁) increases as a function of worsening lung injury (Mori et al., 2018). In this sense, the present study complements our previous model analysis since together they provide a complete analytic description of how a deep inflation of the injured lung recruits closed lung units and then how these units derecruit again over time as governed by the values of the three fixed parameters α, s_m, and k, where α and k depend on V_t. The model of the present study accurately predicts the manner in which the temporal evolution of D_rate varies with both time and injurious V_t (Figure 1), dependencies that are embodied in a single parameter, k, thought to largely reflect how surface tension increases with time (Figure 2). Surface tension elevation is a consequence of surfactant dysfunction, itself caused either by direct damage to the surfactant or by the presence of material that has leaked from the vasculature into the airspaces of the lung. H can also increase as a result of accumulated fluid that floods alveoli. This model thus links developments in physical injury of the parenchyma tissue to the dynamics of lung function derangement via the physiologic mechanism of airspace derecruitment.

Analytical models generally do not have the same flexibility as numerical models to describe the complex details of experimental data. On the other hand, when the model equations are based on putative underlying mechanisms, as in the present case, analytical models may offer insight of a general nature that more empirical nature of numerical models cannot. In this regard, the model in the present study supports the suitability of a power-law function as describing the velocity distribution along virtual trajectories (Eq. 4). The parameter s_m sets the minimum closing velocity, which corresponds to those units that close so slowly they remain open throughout the entire derecruitability test. It is possible that these units would eventually close if the derecruitability tested lasted long enough, but they make no contribution to the data when the tests last only 3 min. Interestingly, the exponent of the power law was found to be −2 purely on the basis of the empirically observed linear increase in lung stiffness following a recruitment maneuver at zero PEEP. In other words, this exponent value comes from the qualitative behavior of the data rather than the fit of a particular model. This is significant because power laws have been invoked in a previous model of lung recruitment during lung inflation to account for avalanches of sequential opening events along the airway tree that are triggered when critical opening pressures are surpassed (Suki et al., 1994, 2000). Moreover, both inter-crackle intervals and terminal air space sizes that are linked by avalanches throughout the airway tree are characterized by power-law distributions having an exponent of −2 (Suki et al., 2000; Alencar et al., 2001). The similarity between the probability distribution for s_c and these other avalanche-based mechanisms in the lung might not be a coincidence, but rather might follow in some way from the fractal branching structure of the airway tree. Thus, although our model consists of a parallel collection of respiratory units that cannot open in the kinds of cascades observed in the airway tree (Suki, 2002), the fractal structure of the lung is embedded in our model through the power-law form of the distribution of s_c (Eq. 4) and thus conforms to previous work by Suki et al. (1994, 2000), Alencar et al. (2001), and Suki (2002).

Having the geometry of the airway tree be implicit in our model in the functional form for s_c could also explain the findings of Massa et al. (2008) that f_so(s₀) and f_sc(s_c) are unchanged following lung injury. That is, although lung injury would be expected to damage the alveolar parenchyma and the airway epithelium, there is no reason to suppose it would alter the branching structure of the airway tree, which would leave the functional form for s_c preserved. On the other hand, we do expect that injury would increase the mean values of the distributions for P_c and P_o (Eqs 2 and 3, respectively) on the basis of increases in surface tension at the air-liquid interface within the lungs as a result of blood-gas barrier perforation and the subsequent airspace accumulation of fluid and protein (Massa et al., 2008; Smith et al., 2013). This conveniently allows us to account for the progression of lung injury largely in terms of the evolution of the value of a single parameter, k in Eq. 9 (Figures 2, 3).

We previously reported an apparent master relationship between and D_rate and H₁ (Smith et al., 2013). The model developed in the present study, however, suggests this is not quite the case because rate of change of D_rate with injury (Eq. 10) is predicted to depend on k, which depends on V_t (Figure 2B), even though the initial value of D_rate (Eq. 10 of the Supplementary Material) does not (Figure 2A). On the other hand, the dependence of dD_rate/dH₀ on V_t is rather subtle (Figure 4), and is easily masked by noise in the data, so the notion of a single master relationship may still be a useful approximation. Also, both H₁ (Mori et al., 2018) and D_rate (Figure 4) exhibit exponential growth over the time-scale of our experiment (up to 4 h), but they are unlikely to continue to increase indefinitely beyond this point because this would correspond to derecruitment of the entire lung. Even if this was a possibility, however, forcing a sizeable V_t into an ever shrinking lung would eventually lead to pressures that would cause rupture of the tissues. In any case, it appears that neither eventually would ever be realized because the mice died when about 75% of the lung became derecruited, at which point H₁ and D_rate were still increasing exponentially.

Another feature of the relationship between D_rate and H₁ is the early transient seen in Figure 4, which shows that D_rate to begin to increase early in the development of lung injury before H₁ changes noticeably. We have previously postulated (Smith et al., 2013) that increases in D_rate are due to increases in surface tension that give rise to derecruitment dynamics manifesting over a timescale of minutes, while increases in H₁ reflect flooded alveoli that have become permanently derecruited. This being the case, the relationship between D_rate and H₁ in Figure 4 suggests that surfactant deactivation by plasma-derived fluid and proteins starts to occur as soon as this material arrives in the airspaces whereas a certain volume of fluid as to accumulate in the airspaces before it can cause noticeable alveolar flooding. It also suggests that D_rate is more directly related to the “rich-get-richer” hypothesis of VILI that we have previously advanced (Hamlington et al., 2018a; Mori et al., 2018) and thus may be the more sensitive marker of the onset of lung injury. On the other hand, even though the model errors for D_rate do not appear to depend on H₁ (Figure 5), at large t, D_rate, and H₁ become linearly dependent on each other (Eq. 10) suggesting that both variables ultimately embody the same pathophysiological information, although the fits of the present study are not as good at this end of the data range compared to those of our previous study (Smith et al., 2013) which curve over and follow the data better. Regardless, both D_rate and H₁ are increased by fluid accumulation in the airspaces, so some kind of monotonic relationship between the two parameters is to be expected.

Our model has a number of limitations. For example, although we have characterized derecruitability of the lung in terms of only two parameters, H₁ and D_rate, the increase in H with time during a derecruitability test is only approximately linear in its early stages; eventually, H begins to plateau (Allen et al., 2002; Allen and Bates, 2004; Massa et al., 2008). Such a plateau is inevitable as only those lung units for which P < P_c will eventually close as τ→∞. D_rate is thus an empirical characterization of only part of a more complex relationship. We also neglected intra-breath variations in the pressures applied to the lungs by mechanical ventilation, as well as the way they depend on the degree of lung injury. We assumed a constant pressure for derecruitment in order to be able to solve the model equations analytically (Eq. 7) and so gain the benefits of insight that a closed-form solution provides. In reality of course, these pressures are not constant because mean airway pressure rises for a given V_t as the lung becomes progressively more consolidated, so more quantitatively accurate predictions would require the numerical solution of a computational model in which airway pressures are allowed to vary in a realistic manner. This may partly explain why the ellipse areas of the joint parameter confidence regions in Figure 3 tend to be greater for larger V_t (i.e., greater VILI); we would generally expect a simple analytical model to fit better to data from normal or mildly injured lungs compared to those that are severely injured. We also waived details concerning alveolar micromechanics such as interactions between neighboring respiratory units, although we recognize the potentially important role of such interactions in lung injury (Knudsen and Ochs, 2018). Finally, we have assumed that the dynamics of derecruitability in the lung are accurate indicators of the degree of lung injury, but they are not direct measures of such. Nevertheless, numerous prior studies from our laboratory have established a tight link between derecruitability and direct measures of damage to the blood-gas barrier (Allen and Bates, 2004; Allen et al., 2005; Smith et al., 2017; Hamlington et al., 2018a,b).

In summary, we have developed an analytical version of our previously developed computational model of the dynamics of lung derecruitment following a recruitment maneuver. We found that the distribution of velocity constants defining rates of movement along the virtual trajectories associated with the pulmonary airways follows a power-law with the same exponent of −2 as found in other investigations into the dynamics of lung recruitment and derecruitment. We speculate that this may reflect the fractal-like structure of the airway tree throughout which recruitment and derecruitment may take place as cascades of open and closing events. Our study also suggests that the rate of change in lung stiffness represented by the parameter D_rate, to the extent that this remains constant during a 3 min derecruitability test, is a particularly sensitive marker of the degree of lung injury.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Ethics Statement

The animal study was reviewed and approved by Institutional Animal Care and Use Committee of the University of Vermont.

Author Contributions

VM developed the model analysis, coded the simulations, and drafted the manuscript. BJS, BS, and JB contributed to the development of the model concept, checked the mathematics, and edited and approved the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was funded by NIH grants R01 HL-124052, R00 HL-128944, and R01 HL-142702 and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brazil (CAPES) – Finance Code 88881.135413/2016-01. These funding agencies contributed to the salary support of the authors of this study. R01 HL-142702 will provide the publication fees for this manuscript.

Conflict of Interest

BJS and JB filed the following patent that is potentially related to the subject matter of the manuscript: “Variable ventilation as a diagnostic tool for assessing lung mechanical function” JB and BJS. PCT Application WO2015127377 A1, Filed on February 23, 2014 (C538).

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2020.542744/full#supplementary-material

References

Albert, S. P., Dirocco, J., Allen, G. B., Bates, J. H., Lafollette, R., Kubiak, B. D., et al. (2009). The role of time and pressure on alveolar recruitment. J. Appl. Physiol. 106, 757–765. doi: 10.1152/japplphysiol.90735.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Alencar, A. M., Buldyrev, S. V., Majumdar, A., Stanley, H. E., and Suki, B. (2001). Avalanche dynamics of crackle sound in the lung. Phys. Rev. Lett. 87: 088101.

Google Scholar

Allen, G., and Bates, J. H. (2004). Dynamic mechanical consequences of deep inflation in mice depend on type and degree of lung injury. J. Appl. Physiol. 96, 293–300. doi: 10.1152/japplphysiol.00270.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, G., Lundblad, L. K., Parsons, P., and Bates, J. H. (2002). Transient mechanical benefits of a deep inflation in the injured mouse lung. J. Appl. Physiol. 93, 1709–1715. doi: 10.1152/japplphysiol.00473.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, G. B., Pavone, L. A., Dirocco, J. D., Bates, J. H., and Nieman, G. F. (2005). Pulmonary impedance and alveolar instability during injurious ventilation in rats. J. Appl. Physiol. 99, 723–730. doi: 10.1152/japplphysiol.01339.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Bates, J. H., and Irvin, C. G. (2002). Time dependence of recruitment and derecruitment in the lung: a theoretical model. J. Appl. Physiol. 93, 705–713. doi: 10.1152/japplphysiol.01274.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Dos Santos, C., and Slutsky, A. (2000). Invited review: mechanisms of ventilator-induced lung injury: a perspective. J. Appl. Physiol. 89, 1645–1655. doi: 10.1152/jappl.2000.89.4.1645

PubMed Abstract | CrossRef Full Text | Google Scholar

Dreyfuss, D., and Hubmayr, R. (2016). What the concept of VILI has taught us about ARDS management. Intens. Care Med. 42, 811–813. doi: 10.1007/s00134-016-4287-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Dreyfuss, D., and Saumon, G. (1998). Ventilator-induced lung injury: lessons from experimental studies. Am. J. Respir. Crit. Care Med. 157, 294–323. doi: 10.1164/ajrccm.157.1.9604014

PubMed Abstract | CrossRef Full Text | Google Scholar

Gattinoni, L., Carlesso, E., Cadringher, P., Valenza, F., Vagginelli, F., and Chiumello, D. (2003). Physical and biological triggers of ventilator-induced lung injury and its prevention. Eur. Respir. J. 22, 15s–25s.

Google Scholar

Gattinoni, L., Marini, J. J., Pesenti, A., Quintel, M., Mancebo, J., and Brochard, L. (2016). The “baby lung” became an adult. Intens. Care Med. 42, 663–673. doi: 10.1007/s00134-015-4200-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Gattinoni, L., and Pesenti, A. (2006). “The concept of baby lung,” in Applied Physiology in Intensive Care Medicine, (Berlin: Springer), 303–311. doi: 10.1007/3-540-37363-2_44

CrossRef Full Text | Google Scholar

Hamlington, K. L., Bates, J. H., Roy, G. S., Julianelle, A. J., Charlebois, C., Suki, B., et al. (2018a). Alveolar leak develops by a rich-get-richer process in ventilator-induced lung injury. PLoS One 13:e0193934. doi: 10.1371/journal.pone.0193934

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamlington, K. L., Smith, B. J., Dunn, C. M., Charlebois, C. M., Roy, G. S., and Bates, J. H. T. (2018b). Linking lung function to structural damage of alveolar epithelium in ventilator-induced lung injury. Respir. Physiol. Neurobiol. 255, 22–29. doi: 10.1016/j.resp.2018.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Knudsen, L., and Ochs, M. (2018). The micromechanics of lung alveoli: structure and function of surfactant and tissue components. J. Histochem. Cell Biol. 150, 661–676. doi: 10.1007/s00418-018-1747-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Massa, C. B., Allen, G. B., and Bates, J. H. (2008). Modeling the dynamics of recruitment and derecruitment in mice with acute lung injury. J. Appl. Physiol. 105, 1813–1821. doi: 10.1152/japplphysiol.90806.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Mori, V., Smith, B. J., Suki, B., and Bates, J. H. (2018). Linking physiological biomarkers of ventilator-induced lung injury to a rich-get-richer mechanism of injury progression. Ann. Biomed. Eng. 47, 638–645. doi: 10.1007/s10439-018-02165-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Parker, J. C., Hernandez, L. A., and Peevy, K. J. (1993). Mechanisms of ventilator-induced lung injury. Crit. Care Med. 21, 131–143. doi: 10.3109/9781420019261-9

CrossRef Full Text | Google Scholar

Robertson, H. T., Altemeier, W. A., and Glenny, R. W. (2000). Physiological implications of the fractal distribution of ventilation and perfusion in the lung. Ann. Biomed. Eng. 28, 1028–1031. doi: 10.1114/1.1310215

CrossRef Full Text | Google Scholar

Schuessler, T. F., and Bates, J. H. (1995). A computer-controlled research ventilator for small animals: design and evaluation. IEEE Trans. Biomed. Eng. 42, 860–866. doi: 10.1109/10.412653

CrossRef Full Text | Google Scholar

Seah, A. S., Grant, K. A., Aliyeva, M., Allen, G. B., and Bates, J. H. (2011). Quantifying the roles of tidal volume and PEEP in the pathogenesis of ventilator-induced lung injury. Ann. Biomed. Eng. 39, 1505–1516. doi: 10.1007/s10439-010-0237-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Slutsky, A. S., and Ranieri, V. M. (2000). Mechanical ventilation: lessons from the ARDSNet trial. Respirat. Res. 1:2.

Google Scholar

Slutsky, A. S., and Ranieri, V. M. (2013). Ventilator-induced lung injury. New Engl. J. Med. 369, 2126–2136.

Google Scholar

Smith, B. J., Bartolak-Suki, E., Suki, B., Roy, G. S., Hamlington, K. L., Charlebois, C. M., et al. (2017). Linking ventilator injury-induced leak across the blood-gas barrier to derangements in murine lung function. Front. Physiol. 8:466. doi: 10.3389/fphys.2017.00466

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, B. J., and Bates, J. H. (2013). Assessing the progression of ventilator-induced lung injury in Mice. IEEE Trans. Biomed. Eng. 60, 3449–3457. doi: 10.1109/tbme.2013.2267151

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, B. J., Grant, K. A., and Bates, J. H. (2013). Linking the development of ventilator-induced injury to mechanical function in the lung. Ann. Biomed. Eng. 41, 527–536. doi: 10.1007/s10439-012-0693-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, B. J., Roy, G. S., Cleveland, A., Mattson, C., Okamura, K., Charlebois, C. M., et al. (2020). Three alveolar phenotypes govern lung function in murine ventilator-induced lung injury. Front. Physiol. 11:660. doi: 10.3389/fphys.2020.00660

CrossRef Full Text | Google Scholar

Suki, B. (2002). Fluctuations and power laws in pulmonary physiology. Am. J. Respirat. Crit. Care Med. 166, 133–137. doi: 10.1164/rccm.200202-152pp

PubMed Abstract | CrossRef Full Text | Google Scholar

Suki, B., Alencar, A. M., Tolnai, J., Asztalos, T., Peták, F., Sujeer, M. K., et al. (2000). Size distribution of recruited alveolar volumes in airway reopening. J. Appl. Physiol. 89, 2030–2040. doi: 10.1152/jappl.2000.89.5.2030

PubMed Abstract | CrossRef Full Text | Google Scholar

Suki, B., Barabási, A.-L., Hantos, Z., Peták, F., and Stanley, H. E. (1994). Avalanches and power-law behaviour in lung inflation. Nature 368:615. doi: 10.1038/368615a0

PubMed Abstract | CrossRef Full Text | Google Scholar

West, B. J. (1990). Physiology in fractal dimensions: error tolerance. Ann. Biomed. Eng. 18, 135–149. doi: 10.1007/bf02368426

PubMed Abstract | CrossRef Full Text | Google Scholar

West, G. B., Brown, J. H., and Enquist, B. J. (2000). “The origin of universal scaling laws in biology,” in Scaling in Biology, eds J. H. Brown and G. B. West (Oxford: Oxford University Press), 87–112.

Google Scholar