On the Stability of Lung Parenchymal Lesions with Applications to Early Pneumothorax Diagnosis

Abstract

Spontaneous pneumothorax, a prevalent medical challenge in most trauma cases, is a form of sudden lung collapse closely associated with risk factors such as lung cancer and emphysema. Our work seeks to explore and quantify the currently unknown pathological factors underlying lesion rupture in pneumothorax through biomechanical modeling. We hypothesized that lesion instability is closely associated with elastodynamic strain of the pleural membrane from pulsatile air flow and collagen-elastin dynamics. Based on the principles of continuum mechanics and fluid-structure interaction, our proposed model coupled isotropic tissue deformation with pressure from pulsatile air motion and the pleural fluid. Next, we derived mathematical instability criteria for our ordinary differential equation system and then translated these mathematical instabilities to physically relevant structural instabilities via the incorporation of a finite energy limiter. The introduction of novel biomechanical descriptions for collagen-elastin dynamics allowed us to demonstrate that changes in the protein structure can lead to a transition from stable to unstable domains in the material parameter space for a general lesion. This result allowed us to create a novel streamlined algorithm for detecting material instabilities in transient lung CT scan data via analyzing deformations in a local tissue boundary.

1. Introduction

Spontaneous pneumothorax is a form of sudden lung collapse closely associated with risk factors such as lung cancer and emphysema. The disorder represents a pressing clinical challenge, affecting between 20 and 40% of patients with major trauma [1]. The onset of pneumothorax has been attributed to the rupture of quasi spherical parenchymal lesions or air blebs of the visceral pleural membrane [2]. The rupture of these lesions in turn releases air into the pleural cavity, and the consequent pressure buildup leads to the collapse of the lung [3]. Because pneumothorax has a relatively large recurrence rate of approximately 54% in the first four years after surgery, wholly understanding and quantifying the mechanisms behind lesion rupture are important [4].

The etiology and trigger mechanisms behind lesion rupture, however, still remain elusive [2, 5]. We hypothesize that lesion rupture and instability are closely associated with elastodynamic strain of the visceral pleural membrane from pulsatile air flow and changes in the constitutive protein composition of tissue. In order to assess the validity of this hypothesis, we constructed a biomechanical model based on the principles of finite strain continuum mechanics and fluid-structure interaction. Our proposed model closely aligns with the behavior of true biological soft tissue through the coupling of membrane energy limitations and growth and remodeling driven by collagen-elastin dynamics. Through the use of an exhaustive ordinary differential equation stability analysis, we isolated several instability regions in the material parameter space of a general lesion.

Based on these mathematical results, we developed an algorithm to rapidly assay clinical lung CT scan data for structural instabilities. The algorithm is based on cataloging and processing local deformations in a tissue boundary isolated via Sobel Edge Detection. The local deformation of the boundary is then fit to our biomechanical model and registered as either stable or unstable based on whether or not the interpolated material parameters lie in the instability region. Figure 1 shows a general outline of our novel early diagnostic method. In the future, the algorithm may be real time coupled with either ultrasound or CT scan data to quantify patient risk of pneumothorax earlier and improve the current diagnostic benchmark of only 75% lesion sensibility [6]. Unlike competitive diagnostic procedures for pneumothorax such as infrared thermography [1], our proposed method relies on detecting pathological hallmarks in tissue behavior before acute lung collapse actually occurs.

Figure 1 

Outline of novel early diagnosis method. An outline of our novel diagnosis method. Based on the observation of certain symptoms, changes in the lung boundary are catalogued and material parameters are derived. The differential equation system is locally applied, and then a stability calculation is performed to determine the appropriate course of action.

2. Methodology

Our study focused on the construction of a theoretically based biomechanical model for lung parenchymal lesion rupture and then the implementation of this model towards a streamlined algorithm for early pneumothorax diagnosis. We began with the theoretical development of a model for simple isotropic tissue deformation and then conducted several initial numerical investigations. Next, we derived mathematical instabilities for our system and then translated these mathematical instabilities to physically relevant structural instabilities via the incorporation of an energy limiter. We then incorporated collagen-elastin dynamics (dynamic alteration of the constitutive proteins in the tissue) to demonstrate how changes in tissue structure may lead to lesion rupture. Lastly, we developed an algorithm for determining lesion stability from CT scan data by fitting our biomechanical model onto tissue deformation data.

2.1. Mathematical Formulation

We modeled the parenchymal lesions with a spherical membrane geometry with thickness significantly smaller than the radius. Several parameters were needed in order to fully characterize the lesion geometry. The nondimensional stretch ratio 𝜆(𝑡) was defined as 𝑟(𝑡)/𝑅, where 𝑟(𝑡) represents the deformed radius varying with time and 𝑅 is the original undeformed radius. Assuming membrane incompressibility, the constant volume condition 𝑉=4𝜋𝑟²(𝑡)ℎ(𝑡) implies that the deformed thickness is ℎ(𝑡)=𝐻/𝜆(𝑡)², where 𝐻 represents the original undeformed thickness. Figure 2 depicts the geometry of these quasi-spherical lesions, including 𝑟(𝑡) and ℎ(𝑡).

Figure 2 

Schema of spherical parenchymal lesion. A schematic drawing of the spherical membrane geometry used as the basis for studying the growth of axisymmetric lung parenchymal lesions.

2.1.1. Constitutive Parenchymal Wall Model

Based on the findings of Tai and Lee [7], who concluded that lung tissue exhibits less than 10% anisotropy of the mean deformation, we pursued an isotropic model for the lung parenchymal wall. We had based a collagen-only isotropic pseudostrain-energy function for the wall of the parenchymal lesion on the work of Denny and Schroter [8]. Previous works by Denny and Schroter [9, 10] had also considered collagen and elastin using separate models [9], which were enhanced with viscoelastic contributions [10]. We adapted the Denny-Schroter constitutive relationship for our 1D formulation through multiplication by the deformed thickness ℎ: 𝐻𝑤=𝜆(𝑡)²𝑐₁𝑒log1−^𝐸11−1𝑐₂+𝑐₃𝐸₁₁, (1) where 𝑐₁ and 𝑐₃ are constants with dimensions N·m⁻², 𝑐₂ is dimensionless, and 𝐸₁₁=(𝜆(𝑡)²−1)/2 (where 𝐄 is the Green-Lagrangian strain tensor). Note that for simplicity, our model does not account for large stiffnesses for large strains. However, this is something that we hope to incorporate in future work. In order to facilitate the stability analysis and further calculations, this pseudostrain-energy function was approximated by a Taylor series centered about 𝐸₁₁=0, yielding −𝑐𝑤=¹𝐻𝑐₂+𝑐₃𝐻𝐸₁₁+−𝑐₁𝐻2𝑐²₂+3𝑐₁𝐻2𝑐₂−2𝑐₃𝐻𝐸²₁₁+𝑂𝐸³₁₁. (2) If the 2D deformation gradient 𝐅 is diag(𝜆(𝑡),𝜆(𝑡)), then the Cauchy stress resultant tensor 𝐓 for inner membrane stress can be represented by the following equation [11]: 𝑇_𝛼𝛽=1𝐹det𝐅^𝛼𝛾𝐹_𝛽𝛿𝜕𝑤𝜕𝐸_𝛾𝛿+2𝜇_𝑚𝐻𝜆(𝑡)³𝑑𝜆(𝑡)𝑑𝑡,𝛼,𝛽,𝛾,𝛿=1,2. (3) Thus, the isotropic stress resultant 𝑇 (𝑇₁₁) is given by the following: 𝑐𝑇=−¹𝐻𝑐₂+𝑐₃−𝑐𝐻+¹𝐻2𝑐²₂+3𝑐₁𝐻2𝑐₂−2𝑐₃𝐻𝜆(𝑡)²+−12𝜇_𝑚𝐻𝜆(𝑡)³𝑑𝜆(𝑡).𝑑𝑡 (4)

2.1.2. Model for Breathing Motion

We used a Fourier series developed by Jakuszkin [12] to model the pulsatile nature of breathing. The driven flow value of air through the lungs 𝑄 is 𝑄=⁷_𝑖=1𝐴_𝑖cos2𝜋𝑓_𝑖𝑡+𝜙_𝑖. (5) Table 1 contains the constants 𝐴_𝑖, 𝑓_𝑖, and 𝜙_𝑖 [12].

Table 1 

Lung Fourier series parameter table.

Note that while we have employed a sinusoidal model for simplicity, we can extend the model to account for erratic breathing behavior, which is not considered here. The driven flow value is related to the pulsatile pressure through the relation Δ𝑃=𝑅_𝑏𝑄, where 𝑅_𝑏 was the respiratory resistance of the bronchiole network. Because 𝑅_𝑏 is 1.0 mmHg·s·L⁻¹ [13] and the original pressure was atmospheric pressure (760 mmHg), the final pressure series was as follows: 𝑃_inner=760+⁷_𝑖=1𝐴_𝑖cos2𝜋𝑓_𝑖𝑡+𝜙_𝑖. (6)

2.1.3. Dynamic Equation System

We obtained the dynamic equation for the wall through linear momentum balance to yield 𝜌_𝑚𝑑ℎ𝑅²𝜆(𝑡)𝑑𝑡²(𝑡)−(𝜆(𝑡=𝑃2𝑇))𝜆(𝑡)𝑅, (7) where 𝜌_𝑚 is the density of the membrane and 𝑃(𝑡) is the inner pressure minus the outer pressure (transmural pressure). The Navier-Stokes equations in spherical coordinates were used to model the pressure that the pleural fluid exerted on the lung parenchyma. The resultant pressure was as follows [14]: 𝑃_outer=𝑝_∞+𝜌_𝑓𝑅²𝜆(𝑡)𝑑²𝜆(𝑡)𝑑𝑡²+32𝑑𝜆(𝑡)𝑑𝑡², (8) where 𝑝_∞ is a constant value representing the pressure at a significantly large distance from the lung and 𝜌_𝑓 is the density of the fluid. Note the absence of a viscous term as well as the same convention for the pressure at infinity as in the published work by Shah and Humphrey [14]. Combining the equations for the dynamic wall, the outer pressure, the dynamic radial stress which contributes the factor of 𝜇_𝑓, and viscoelasticity, we arrived at the following differential equation for the air-tissue-pleural fluid system: 𝜌_𝑚𝐻𝑅𝜆(𝑡)²+𝜌_𝑓𝑅²𝜆(𝑡)𝑑²𝜆(𝑡)𝑑𝑡²+32𝜌_𝑓𝑅²𝑑𝜆(𝑡)𝑑𝑡²𝜇+4^𝑓𝜆(𝑡)×𝑑𝜆(𝑡)+(𝜆(𝑡𝑑𝑡2𝑇))𝜆(𝑡)𝑅+4𝜇_𝑚𝐻𝜆(𝑡)⁴𝑅𝑑𝜆(𝑡)𝑑𝑡=𝑃_inner(𝑡)−𝑝_∞, (9) where 𝜇_𝑓 and 𝜇_𝑚 are the dynamic viscosities of the pleural fluid and membrane, respectively.

2.1.4. Nondimensionalization

Nondimensionalization is a process that helps to accommodate multiscale variations in the data. This technique is a common engineering practice which essentially simplifies multiphysics problems with different measured units involved. Because our final governing equations couple three different physical aspects (the air, the membrane, and the pleural fluid), this method was necessitated. Moreover, nondimensionalization recovers the characteristic properties of the stretch ratio and stretch rate. We followed this standard which is commonly used in published works on biomechanics such as Shah and Humphrey [14]. The following nondimensionalized quantities were required: 𝑦₀≡𝜆(𝑡),𝜏≡𝑡𝑐₁𝐻𝜌_𝑚𝑅²𝐻𝜌,𝑏≡^𝑓𝑅𝜌_𝑚𝐻,𝐻𝑑≡𝑅,𝑚_𝑓≡𝜇_𝑓𝜌_𝑚𝑐₁𝐻²𝜇,𝑚≡^𝑚𝜌_𝑚𝑐₁𝐻²,𝐹(𝜏)≡𝑃_inner(𝑡)−𝑝_∞𝑅𝑐₁𝐻𝑇,𝑓≡𝑐₁𝐻. (10) The resulting nondimensionalized equation system was as follows: 𝑑𝑦₀𝑑𝜏=𝑦₁,𝑑𝑦₁=𝐹(𝜏)𝑑𝜏−3𝑏𝑦²₁/2−4𝑚_𝑓𝑦₁/𝑦₀−2𝑓𝑦₀/𝑦₀−4𝑚𝑑𝑦₁/𝑦⁴₀𝑏𝑦₀+𝑦₀⁻². (11)

2.1.5. Parameter Table

Table 2 is a complete collection of the relevant parameter values for the model.

Table 2 

List of parameters.

2.2. Initial Numerical Investigations

We conducted several numerical investigations in order to quantify both system dynamics and stability. The previous dynamical system was numerically evolved in MATLAB through a fourth-order explicit Runge-Kutta scheme. Plots of the stretch ratio 𝜆(𝑡) versus time reveal sustained stable oscillations of the membrane for normal lung parenchymal tissue parameters of the Denny-Schroter model: 𝑐₁=−22.5×10⁵, 𝑐₂=1.26, and 𝑐₃=−7.8×10⁵. Variations in material stiffness proved an effective tool in analyzing the effects of parameter variation on system dynamics. Figure 3 shows that as material stiffness increased in the model (for test values of 𝑐₂=1.00,1.26,and1.52), the amplitude of the oscillations decreased as one would expect intuitively. Figure 4 shows that for the general range of 𝑐₂∈[.5,1.6], amplitude can be seen to monotonically decrease and begin to level off near zero. For the same window of stiffness constants, oscillation frequency reaches a peak near 𝑐₂≈0.78 and then decreases monotonically from that value.

Figure 3 

Stretch factor 𝜆(𝑡) for different stiffness constant (𝑐2). This figure shows three plots of stretch factor with respect to dimensionless time for the stiffness constants 𝑐2=1.00,1.26,and1.52. Variations in the stiffness constant clearly affect the amplitude and frequency of oscillations.

(a)

(b)

Figure 4 

Quantified effects of stiffness constant (𝑐2) on amplitude and frequency of oscillations. (a) The graph shows the effect of variations in stiffness constant on the amplitude of oscillations. The monotonically decreasing nature of the graph suggests that the stiffer the membrane, the smaller the amplitude. (b) The graph shows the effect of variations in stiffness constant on the frequency of oscillations with a maximum at 𝑐2≈0.78.

Figure 5 demonstrates that the system is stable even with minor perturbations in initial conditions for that class of material parameters. This suggests that there may be conditions where the bleb will not rupture. Figure 6 is explained in the preceding paragraph as a relative force proportion graph, which illustrates the dynamic behavior of the different forces including internal pressure, fluid structure, internal membrane, radial stress, and viscoelasticity. The graph suggests that friction is not a significant factor leading to pneumothorax, and this aligns with other published works. The numerical phase plane analysis in Figure 5 of 𝜆(𝑡) versus 𝜆̇(𝑡) further suggests the existence of an orbit about a stable equilibria. The spiral nature of the superimposed vector field, representing the stretch rate 𝜆̇(𝑡) versus the derivative of the stretch rate 𝜆̈(𝑡), suggests an asymptotically stable node at 𝜆_eq≈1.3 for normal lung parenchymal tissue parameters. A relative force proportion graph in Figure 6 reveals the forces of internal pressure and inner membrane stress as the major driving forces behind this oscillation, whereas the viscoelastic and fluid forces effects are negligible. Furthermore, several unstable regions had been isolated through variation of the material parameters such as for 𝑐₁=−22.5×10⁵, 𝑐₂=1.26, and 𝑐₃=−1.39×10⁶, but a rigorous mathematical analysis had been called for to fully determine whether the instabilities were truly physical or just due to numerical effects.

(a)

(b)

Figure 5 

Numerically generated phase planes of the dynamical system. (a) The stable orbit and superimposed vector field demonstrated in the phase plane suggest the presence of an asymptotically stable equilibrium. (b) The 3D plot depicts the evolution of the stable orbit with respect to dimensionless time as further evidence of sustained stable oscillations. Refer to Section 2.4 for more information about physically unstable situations.

(a)

(b)

Figure 6 

Relative force proportion graph. (a) The graph depicts the relative magnitude of the dimensionless forces at play over dimensionless time: pressure (𝐹(𝜏)), radial stress (4𝑚𝑓𝑦1/𝑦0), internal membrane (2𝑓(𝑦0)/𝑦0), fluid structure (3𝑏𝑦21/2), and viscoelastic (4𝑚𝑑𝑦1/𝑦0). It is clear that internal pressure and internal membrane forces are the major driving forces behind the oscillations. (b) The three forces that contribute the least (radial stress, viscoelastic, and fluid structure) are plotted on semilogarithmic axes.

2.3. Stability Results

2.3.1. Determining System Equilibria

In order to simplify the following equilibria identification and stability analysis, we made the system autonomous by approximating 𝐹(𝑡)≈𝐹=𝑅𝑃/𝑐₁𝐻. The system is at an equilibrium for any point (𝜆(𝑡),𝜆̇(𝑡)) such that 𝑦̇₀=0 and 𝑦̇₁=0. This implies that 𝜆̇(𝑡)=0 for any equilibrium point. The equilibrium stretch ratio 𝜆_eq may be found by solving 𝐹−(𝑓(𝜆(𝑡))/𝜆(𝑡))=0 as follows: 𝜆_eq=−𝛼±−4𝛽𝛾+𝛼²,2𝛾 (12) where 𝛼=𝑐₁𝑐²₂𝐹, 𝛽=−𝑐₁+5𝑐₁𝑐₂−6𝑐²₂𝑐₃, and 𝛾=𝑐₁−3𝑐₁𝑐₂+4𝑐²₂𝑐₃. For physically meaningful values of 𝜆(𝑡), the positive equilibrium point is the one of significance. For the normal lung parenchymal tissue parameters, 𝜆⁺_eq=1.3405 as suggested by the initial numerical investigation.

2.3.2. Characterization of Unstable Nodes

We conducted a general ordinary differential equation (ODE) stability analysis on the physical system in order to evaluate a lung parenchymal lesion with arbitrary material parameters for physical instabilities. Given a system with Jacobian matrix 𝐽, ODE instability is prescribed for det𝐽<0ortr𝐽>0. For the previous dynamical system in 𝑦₀ and 𝑦₁, the Jacobian matrix about the positive equilibrium point reduces to𝐽=^{012𝐹𝛾−4𝛽𝛾+𝛼2−𝛼−4𝛽𝛾+𝛼2}_{𝛼−𝛼+−4𝛽𝛾+𝛼224𝛾2/−𝛼+−4𝛽𝛾+𝛼22+𝑏−𝛼+−4𝛽𝛾+𝛼2/2𝛾}^2𝑚_𝑓^{𝛼+−4𝛽𝛾+𝛼2−/𝛽64𝛾4𝑑𝑚/−𝛼+−4𝛽𝛾+𝛼24}_{4𝛾2/−𝛼+−4𝛽𝛾+𝛼22+𝑏−𝛼+−4𝛽𝛾+𝛼2/2𝛾}. (13)

The constraint of physically relevant material parameters imposes the inequalities 𝑐₁<0,  𝑐₂>0, and 𝑐₃<0, whereas the constraint of defined real elements of the Jacobian matrix imposes 𝛽≠0,  𝛾≠0, and 𝛼²≥4𝛽𝛾. Further note that from the definition of 𝐹=𝐴𝑃/𝑐₁𝐻, 𝛼 is always greater than zero for physically relevant material constants (𝐴,𝐻,𝑃>0).

Theorem 1. For any values of 𝛾 and 𝛽 such that 𝛾<0,  𝛽>0, and 𝛾<𝑏^1/3𝛼−𝑏^2/3𝛽, the system is unstable because det𝐽<0.

The expression det𝐽=−𝜕𝑦̇₁/𝜕𝑦₀ is less than zero for all 𝛽,𝛾 such that sgn(𝜕𝑦̇₁/𝜕𝑦₀)>0. The sign of this expression is in turn governed by the following three subexpressions. The subexpression of (−𝛼+−4𝛽𝛾+𝛼²)² is consciously excluded from the sign analysis because it is invariantly positive.

Lemma 2. The sign of 𝛽 and 𝛾 cannot both be positive for physically relevant parameters.
The lemma is proven by contradiction. By definition, 𝛽>0 implies −𝑐₁+5𝑐₁𝑐₂−6𝑐²₂𝑐₃>0. Isolation of 𝑐₃ reduces the inequality to 𝑐₃<𝑐₁(5𝑐₂−1)/6𝑐²₂. Likewise, 𝛾=𝑐₁−3𝑐₁𝑐₂+4𝑐²₂𝑐₃>0 reduces to 𝑐₃>𝑐₁(3𝑐₂−1)/4𝑐²₂ and places bounds on 𝑐₃. However, the inequality 𝑐₁(3𝑐₂−1)/4𝑐²₂<𝑐₁(5𝑐₂−1)/6𝑐²₂ reduces to 0<𝑐₁(𝑐₂+1), which cannot be true from the restrictions of 𝑐₁<0 and 𝑐₂>0.

Lemma 3. The expression −4𝛽𝛾+𝛼²−𝛼−4𝛽𝛾+𝛼² (14) is greater than zero for any value of 𝛽,𝛾 such that either 𝛽<0 and 𝛾>0 or 𝛽>0 and 𝛾<0. Similarly, the expression is less than zero for any value of 𝛽,𝛾 such that either 𝛽<0 and 𝛾<0 or 𝛽>0 and 𝛾>0.
The inequality −4𝛽𝛾+𝛼²−𝛼𝛼²−4𝛽𝛾>0 reduces to 𝛼²−4𝛽𝛾>𝛼𝛼²−4𝛽𝛾. Because both sides of the inequality are positive, the inequality may be squared and the terms brought to one side to give (𝛼²−4𝛽𝛾)(−4𝛽𝛾)>0. Because only values of 𝛼,𝛽,   and 𝛾 such that 𝛼²−4𝛽𝛾>0 are considered, the expression is true 𝛽𝛾<0. Thus, −4𝛽𝛾+𝛼²−𝛼−4𝛽𝛾+𝛼²>0 for either 𝛽<0 and 𝛾>0 or 𝛽>0 and 𝛾<0. Likewise, −4𝛽𝛾+𝛼²−𝛼−4𝛽𝛾+𝛼²<0 for 𝛽<0 and 𝛾<0 or 𝛽>0 and 𝛾>0.

Lemma 4. For physically relevant restriction of 𝑏>0, the expression 4𝛾²−𝛼+−4𝛽𝛾+𝛼²²+𝑏−𝛼+−4𝛽𝛾+𝛼²2𝛾 (15) is greater than zero for any value of 𝛽,𝛾 such that any one of the following three conditions are met: 𝛾>0 and 𝛽<0; 𝛾<0 and 𝛽<0; 𝛾<𝑏^1/3𝛼−𝑏^2/3𝛽,  𝛾<0, and 𝛽>0.
For either the conditions of 𝛾>0 and 𝛽<0 or 𝛾<0 and 𝛽<0, the expression 𝑏(−𝛼+−4𝛽𝛾+𝛼²)/2𝛾 is always greater than zero. Hence, because 4𝛾²/(−𝛼+−4𝛽𝛾+𝛼²)² is invariantly positive, the entire expression (15) will be greater than zero for these conditions. Given 𝛾,𝛽 such that 𝛾<0 and 𝛽>0, the inequality reduces to 8𝛾³+𝑏(−𝛼+−4𝛽𝛾+𝛼²)<0 or 𝛾<−(𝑏^1/3/2)(−𝛼+−4𝛽𝛾+𝛼²). This in turn reduces to 4𝛾/𝑏^2/3−4𝛼/𝑏^1/3+4𝛽<0 or 𝛾<𝑏^1/3𝛼−𝑏^2/3𝛽.

Lemma 5. For physically relevant restriction of 𝑏>0, the expression (15) is less than zero for any value of 𝛽,𝛾 such that 𝛾<0,  𝛽>0, and 𝛾>𝑏^1/3𝛼−𝑏^2/3𝛽.
Either the conditions of 𝛾>0 and 𝛽<0 or 𝛾<0 and 𝛽<0 make the statement unconditionally false, leaving only the case of 𝛾<0,  𝛽>0. The fact that 𝛾>𝑏^1/3𝛼−𝑏^2/3𝛽 is proven likewise as in Lemma 4.

Lemma 6. The expression 2𝐹𝛾𝛼 (16) is greater than zero for all 𝛾<0 and less than zero for all 𝛾>0. By definition, 𝐹<0 and 𝛼>0, and thus the result follows.
The expression det𝐽<0 is true such that the product of expressions (14), (15), and (16) is positive. Upon coupling of the previous inequality systems, the only scenario that does not lead to contradictions is for (14), (15), and (16) all positive. This case implies that for 𝛾<0,  𝛽>0, and 𝛾<𝑏^1/3𝛼−𝑏^2/3𝛽 the system is unstable because det𝐽<0.

Theorem 7. For any values of 𝛾 and 𝛽 such that 𝛾<0, 𝛽>0, 𝛾>𝑏^1/3𝛼−𝑏^2/3𝛽, and 𝛽>𝛼𝛾²(𝑑𝑚/𝑚_𝑓)^1/3+𝛾(𝑑𝑚/𝑚_𝑓)^2/3, the system is unstable because tr𝐽>0.

The expression tr𝐽=−𝜕𝑦̇₁/𝜕𝑦₁ is less than zero for all 𝛽,𝛾 such that sgn(𝜕𝑦̇₁/𝜕𝑦₁)>0. The sign of this expression is in turn governed by subexpressions (15) and (17).

Lemma 8. For physically relevant parameters (𝑚_𝑓,𝑚,𝑑>0), 2𝑚_𝑓𝛼+−4𝛽𝛾+𝛼²𝛽−64𝛾⁴𝑑𝑚−𝛼+−4𝛽𝛾+𝛼²⁴ (17) is less than zero for all 𝛽<0 or for 𝛽>0,  𝛾<0, and 𝛽>𝛼𝛾²(𝑑𝑚/𝑚_𝑓)^1/3+𝛾(𝑑𝑚/𝑚_𝑓)^2/3.
For all 𝛽<0, expression (17) is invariantly negative, and thus, the statements are unconditionally true. This leaves only the case of 𝛽>0 and 𝛾<0. Manipulation of the expression yields 2𝑚_𝑓𝛽𝛼+𝛼²−4𝛽𝛾−𝛼+𝛼²−4𝛽𝛾⁴<64𝛾⁴𝑑𝑚 (18) or (−𝛼+𝛼²−4𝛽𝛾)³>−8𝛾³𝑑𝑚/𝑚_𝑓. This in turn yields 𝛼²−4𝛽𝛾>(𝛼−2𝛾(𝑑𝑚/𝑚_𝑓)^1/3)², which upon isolation of 𝛽 gives the expression 𝛽>𝛼𝛾²(𝑑𝑚/𝑚_𝑓)^1/3+𝛾(𝑑𝑚/𝑚_𝑓)^2/3.

The expression tr𝐽>0 is true such that the product of expressions (15) and (17) is positive. For the signs of (15) and (17) to both be negative, Lemma 5 and Lemma 8 give 𝛾<0, 𝛽>0, 𝛾>𝑏^1/3𝛼−𝑏^2/3𝛽, and 𝛽>𝛼𝛾²(𝑑𝑚/𝑚_𝑓)^1/3+𝛾(𝑑𝑚/𝑚_𝑓)^2/3. The case of expression (15) and (17) both greater than zero is consciously omitted from the proof because the analysis leads to an instability region already enclosed by Theorem 1. Figure 7 shows 2D portraits of the regions obtained from these instability criteria for a certain selection of material parameters 𝑐₁, 𝑐₂, and 𝑐₃.

(a)

(b)

(c)

Figure 7 

2D instability region portraits. The images show the regions of instability (blue) for each selection of two material parameters. The value of the third parameter in each graph is assumed to be the value from Denny and Schroter [8].

2.4. Introducing Energy Limitations of Tissue

In order to accurately capture the biomechanics of lung parenchymal failure, we translated ODE instabilities to physical mechanical failure via the introduction of a finite energy limiter Φ(𝐽𝑚⁻³). Similar to Volokh and Vorp [19], we built a finite energy constitutive model Ψ based on the framework of the hyperelastic model 𝑤 which was constructed as follows: Ψ=𝐻Φ1−exp−𝑤𝐻Φ. (19) Note that as 𝑤→∞, Ψ→𝐻Φ a finite value of energy before rupture in contrast to the hyperelastic model where 𝑤 may grow unbounded as 𝜆(𝑡)→∞. Similarly, in the limiting case 𝐻Φ≪𝑤,  Ψ≈𝑤 as can be seen via Taylor expansion. Reformulating the Cauchy stress resultant tensor 𝐓 given by [14] with the new strain-energy function yields 𝑇^′_𝛼𝛽=1𝐹det𝐅^𝛼𝛾𝐹_𝛽𝛿𝜕Ψ𝜕𝐸_𝛾𝛿+2𝜇_𝑚𝐻𝜆(𝑡)³𝑑𝜆(𝑡)𝑑𝑡,𝛼,𝛽,𝛾,𝛿=1,2. (20) Thus, the new isotropic stress resultant 𝑇^′ (𝑇^′₁₁) is given by the following: 𝑇^′=−𝑐₁𝐻𝑐₂+𝑐₃−𝑐𝐻+¹𝐻2𝑐²₂+3𝑐₁𝐻2𝑐₂−2𝑐₃𝐻𝜆(𝑡)²𝑐−1×exp¹𝐻/𝑐₂−𝑐₃𝐻𝜆(𝑡)²−1/2+𝑐𝐻Φ₁𝐻/2𝑐²₂−3𝑐₁𝐻/2𝑐₂+2𝑐₃𝐻𝜆(𝑡)²−1²/4+𝐻Φ2𝜇_𝑚𝐻𝜆(𝑡)³𝑑𝜆(𝑡).𝑑𝑡 (21)

As demonstrated in Figure 8, numerical simulations with the new isotropic stress resultant demonstrate rupture midexpansion as opposed to expansion to infinity in the original constitutive model.

Figure 8 

Effect of introduction of finite energy limiter. A point (𝑐1,𝑐2,𝑐3) which satisfies the instability of requirement of Theorem 1  (−22.5×105,1.26,and−1.39×106) was used to test the effect of the finite energy limiter. As indicated by the red plot, the finite energy limiter curtails expansion to infinity and establishes a physical rupture point.

2.5. Collagen-Elastin Dynamics

The final component of the proposed model for parenchymal lesion growth and rupture is a means by which stable mechanical configurations may enter unstable domains. The introduction of collagen-elastin dynamics (a dynamic protein framework of the tissue) elucidates mechanisms for how a stable lesion may eventually rupture. We introduced a new constitutive model through the creation of a linear combination of Denny and Schroter’s [8] collagen-only strain energy function and Humphrey and Yin’s [20] elastin-only strain energy function: 𝑤^′=𝐴_𝑐(𝜏𝐻)𝜆(𝑡)²𝑐₁𝑒log1−^𝐸11−1𝑐₂+𝑐₃𝐸₁₁+𝐴_𝑒(𝜏)𝑏𝐻𝜆(𝑡)²𝜆(𝑡)(𝑡),𝑤−log𝜆−1_′≈−𝐴_𝑐(𝜏)𝑐₁𝐻𝑐₂+𝐴_𝑐(𝜏)𝑐₃𝐻𝜆(𝑡)²−12+𝐴_𝑒(𝜏)𝑏𝐻2−𝐴_𝑐(𝜏)𝑐₁𝐻2𝑐²₂+3𝐴_𝑐(𝜏)𝑐₁𝐻2𝑐₂−2𝐴_𝑐(𝜏)𝑐₃𝐻𝜆(𝑡)²−1²4, (22) where 𝐴_𝑐(𝜏) and 𝐴_𝑒(𝜏) are dimensionless quantities representing active collagen and elastin numbers. Similarly as before, we introduced a finite energy limiter Ψ^′ to this constitutive model, and we recalculated the new isotropic strain tensor. Unlike the previous developments of the paper, collagen-elastin dynamics called for a return to the initial momentum balance due to the changing mass of the protein matrix. Given the mass of collagen 𝐶 and elastin 𝐸 defined by 𝜌_𝑐𝐴_𝑐 and 𝜌_𝑒𝐴_𝑒 (where 𝜌_𝑐 and 𝜌_𝑒 are equal to 5.10 kg·m⁻³ [21]), the new momentum balance was as follows: 𝜌_𝑚𝑑ℎ𝑅²𝜆(𝑡)𝑑𝑡²+𝑑𝐴_𝑒𝜌_𝑒+𝐴_𝑐𝜌_𝑐𝑣𝑑𝑡^𝑟(𝑡)−(𝜆(𝑡𝐻=𝑃2𝑇))𝜆(𝑡)𝑅, (23) where 𝑣_𝑟(𝑡) is the velocity at the boundary of an expanding biological membrane given by 𝑅𝑑𝜆(𝑡)/𝑑𝑡. With the nondimensionalized term 𝑚_𝑝=(𝐴/𝜌_𝑚|𝑐₁|)(𝑑𝐶/𝑑𝑡+𝑑𝐸/𝑑𝑡), the first two differential equations of the dynamical system became 𝑑𝑦₀𝑑𝜏=𝑦₁,𝑑𝑦₁=𝐹(𝜏)𝑑𝜏−3𝑏𝑦²₀/2−4𝑚_𝑓𝑦₁/𝑦₀−2𝑓𝑦₀/𝑦₀−4𝑚𝑑𝑦₁/𝑦⁴₀−𝑚_𝑝𝑦₁/𝑦²₀𝑏𝑦₀+𝑦₀⁻². (24) Based on histochemical evidence that had determined that collagen levels increase with increased stress and that elastin levels decrease because of elastolytic processes near rupture [22], we proposed the following differential equations for collagen-elastin dynamics: 𝑑𝐴_𝑐=𝑘𝑑𝑡¹𝐴_𝑐𝜙^′−𝜓^′𝑛,𝑑𝐴_𝑒=𝑑𝑡−𝑘₂𝐴_𝑒𝜙^′−𝜓^′𝑛. (25) Because 𝑘₁ and 𝑘₂ are positive constants, we devised the differential equations such that little change occurs in the protein matrix in scenarios far from rupture (|𝜙^′−𝜓^′|≫0) and larger changes occur as the membrane nears rupture. We numerically solved the series of four differential equations as a whole in MATLAB, revealing escapes from stable orbits for certain rate kinetics of the elastin-collagen dynamics that are depicted in Figure 9.

(a)

(b)

Figure 9 

Collagen-elastin dynamics may lead to escapes from stable orbits. (a) The graph shows 𝜆(𝜏) with the additional collagen-elastin dynamics calibrated with 𝐴𝑐=100,  𝐴𝑒=50,  𝑘1=1×107,  𝑘2=1×103,  and 𝑛=1. The point of rupture occurs at the final peak of oscillation at 𝜏≈16. (b) The phase plane reflects the transition from stability to instability as the trajectory veers off from stable orbits.

2.6. Parameter Estimation from In Vivo Data

In order to build a database for parameter estimation on the proposed model, we collected in vivo lung data in the form of CT scans from the VIA/I-ELCAP Public Lung Image Database, which were about 250 mm × 150 mm in size. The Sobel Edge Detection algorithm provided a means to quantify the deformation of the lung boundary and to extract material parameters from in vivo clinical data for the strain-energy function. This algorithm works by estimating the horizontal and vertical spatial gradients at the interior pixels, 𝐺_𝑥 and 𝐺_𝑦, through application of a matrix mask in order to detect edges. Pixels that meet a certain threshold for 𝐺  (|𝐺_𝑥|+|𝐺_𝑦|) are counted as edges (Aybar, “Sobel Edge Detection”). This method was preferred over methods such as the Canny algorithm due to the coupling of speed and required accuracy for this problem. In order to facilitate the detection of the edges, we recolored the lung images to a 50% black and white scheme with 50% softening. After recoloration and softening, we ran the images through the Sobel algorithm implemented in Python, which returned a list of edge pixels for each image as coordinate pairs. We first transformed these pairs into a new coordinate system with the center of the image as the origin then categorized the pairs into subsections corresponding to sectors defined by angles. We then calculated values for the stretch factor, kinetic energy per square meter ((1/2)𝜌_𝑚𝑣²𝐻), and radial velocity in each radial direction.

Next, we applied nonlinear regression using the Levenberg-Marquadt algorithm on each subsection in order to extract the three material parameters 𝑐₁, 𝑐₂, and 𝑐₃ as well as the corresponding coefficients of determination (𝑅² values). We fit the velocity and kinetic energy data to the Denny-Schroter second-order Taylor series approximation. After we determined these parameters, they were classified as either stable or unstable corresponding to the previously derived stability criteria. Lastly, our program recolored the boundary of the lung parenchyma based on the stability of the corresponding subsections, with green corresponding to stable subsections and red corresponding to unstable subsections. The general schema for our computational algorithm is shown in Figure 10. Thus, our algorithm is streamlined towards interpolating material parameters from deformations in lung CT scan data and identifying specific areas of risk in the lung parenchyma.

(a)

(b)

(c)

Figure 10 

Edge detection and risk analysis. (a) and (b) show a lung before and after deformation, with red pixels indicating edge pixels detected by the Sobel Edge Detection algorithm. (c) shows the results of the nonlinear regression analysis performed on different sectors of the lung. Areas characterized as unstable, and thus at risk for lung collapse, are colored in red, whereas areas characterized as stable are colored in green.

Motion artifacts, while important and are often observed in CT data acquisition, are not considered in this paper. The images are represented using intensity values at various pixel locations. This paper used only one dataset as a proof of concept. Our algorithm looks at small sections of the lung boundary, allowing for the isotropic assumption to be valid. However, it must be noted that this assumption must be removed for more general studies.

3. Conclusion and Future Directions

Our proposed biomechanical model for lung parenchymal lesions demonstrates that elastodynamic strain of the visceral pleural membrane from pulsatile air flow and changes in the constitutive protein composition of tissue are key pathological factors in pneumothorax. Through proving general instability results, our study demonstrated that these factors biject to physically unstable domains in the material parameter space for a general lesion. Specifically, our discovery that certain rate kinetics of collagen-elastin dynamics may lead to unstable mechanical configurations aligns well with clinical studies, which have found that certain connective tissue diseases lead to spontaneous pneumothorax [23, 24]. Future work on the theoretical biomechanics aspect of our model will incorporate COMSOL multiphysics simulation software integrated with finite element mesh processing to solve the problems of nonaxisymmetry and anisotropic lesion growth. Also, it would be interesting to study how variations in the dynamic viscosities affect the behavior of parenchymal lesions. Clinical studies have indicated that gravity is a significant factor on the air accumulated in the pleural cavity. We hope to consider these effects in future work. Quantification of material properties will require efficient inverse parameter estimation studies, which we hope to incorporate in the future. Such studies will help us to develop efficient simulation tools [25, 26]. Other methods for lung deformation registration such as the dictionary learning approach used by Zhang et al. [27] will be considered.

Our algorithm for determining lesion stability from CT scan data relies on fitting our biomechanical model to deformations in a local tissue boundary. Using Sobel Edge Detection and CT scans from a lung image database, we could isolate unstable regions of local lung parenchymal tissue. The algorithm is a robust method for quantifying patient risk of pneumothorax and may be streamlined to work in the future with in vivo lung data collected from either ultrasound or CT scans. Overall, our computational solution for lung parenchymal lesion detection and patient-specific structural instability profiling is a feasible alternative diagnostic strategy and has the potential to surpass current diagnostic benchmarks for pneumothorax.