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

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.


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 fluidstructure 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: 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. 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.

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.

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 2 ( )ℎ( ) implies that the deformed thickness is ℎ( ) = / ( ) 2 , where represents the original undeformed thickness. Figure 2 depicts the geometry of these quasispherical lesions, including ( ) and ℎ( ).

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 pseudostrainenergy 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 ℎ: where 1 and 3 are constants with dimensions N⋅m −2 , 2 is dimensionless, and 11 = ( ( ) 2 − 1)/2 (where E 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 11 = 0, yielding Computational and Mathematical Methods in Medicine If the 2D deformation gradient F is diag( ( ), ( )), then the Cauchy stress resultant tensor T for inner membrane stress can be represented by the following equation [11]: Thus, the isotropic stress resultant ( 11 ) is given by the following:

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 Table 1 contains the constants , , and [12]. 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 −1 [13] and the original pressure was atmospheric pressure (760 mmHg), the final pressure series was as follows:

Dynamic Equation System
. We obtained the dynamic equation for the wall through linear momentum balance to yield 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]: 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: where and are the dynamic viscosities of the pleural fluid and membrane, respectively.

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: The resulting nondimensionalized equation system was as follows: Table. Table 2 is a complete collection of the relevant parameter values for the model.

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: 1 = −22.5 × 10 5 , 2 = 1.26, and 3 = −7.8 × 10 5 . 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 2 = 1.00, 1.26, and 1.52), the amplitude of the oscillations decreased as one would expect intuitively. Figure 4 shows that for the general range of 2 ∈ [.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 2 ≈ 0.78 and then decreases monotonically from that value. 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 ( ) versuṡ ( ) 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 1 = −22.5 × 10 5 , 2 = 1.26, and 3 = −1.39×10 6 , but a rigorous mathematical analysis had been called for to fully determine whether the instabilities were truly physical or just due to numerical effects.

Determining System Equilibria.
In order to simplify the following equilibria identification and stability analysis, we made the system autonomous by approximating ( ) ≈ = / 1 . The system is at an equilibrium for any point where = 1 2 2 , = − 1 + 5 1 2 − 6 2 2 3 , and = 1 − 3 1 2 + 4 2 2 3 . 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.

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 < 0 or tr > 0. ] .
The constraint of physically relevant material parameters imposes the inequalities 1 < 0, 2 > 0, and 3 < 0, whereas the constraint of defined real elements of the Jacobian matrix imposes ̸ = 0, ̸ = 0, and 2 ≥ 4 . Further note that from the definition of = / 1 , is always greater than zero for physically relevant material constants ( , , > 0).

The expression det
= −̇1/ 0 is less than zero for all , such that sgn(̇1/ 0 ) > 0. The sign of this expression is in turn governed by the following three Computational and Mathematical Methods in Medicine 7 subexpressions. The subexpression of (− + √−4 + 2 ) 2 is consciously excluded from the sign analysis because it is invariantly positive.

Lemma 3. The expression
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.

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.
The expression tr = −̇1/ 1 is less than zero for all , such that sgn(̇1/ 1 ) > 0. The sign of this expression is in turn governed by subexpressions (15) and (17). 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 > 2 ( / ) 1/3 + ( / ) 2/3 . The case of expression (15) and (17)  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 1 , 2 , and 3 .

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 Φ( −3 ). 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: 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 T given by [14] with the new strain-energy function yields Thus, the new isotropic stress resultant ( 11 ) is given by the following: 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.

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: 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 −3 [21]), the new momentum balance was as follows: where V ( ) is the velocity at the boundary of an expanding biological membrane given by ( )/ . With the nondimensionalized term = ( /√ | 1 |)( / + / ), the first two differential equations of the dynamical system became 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: Because 1 and 2 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. 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) V 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 1 , 2 , and 3 as well as the corresponding coefficients of determination ( 2 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.
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.

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 patientspecific structural instability profiling is a feasible alternative diagnostic strategy and has the potential to surpass current diagnostic benchmarks for pneumothorax.