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

Authors: Archis R. Bhandarkar; Rohan Banerjee; Padmanabhan Seshaiyer
Journal: Computational and Mathematical Methods in Medicine
 (2013)

---

## 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.


---

## Body

## 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. Figure1 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λ(t) was defined as r(t)/R, where r(t) represents the deformed radius varying with time and R is the original undeformed radius. Assuming membrane incompressibility, the constant volume condition V=4πr2(t)h(t) implies that the deformed thickness is h(t)=H/λ(t)2, where H represents the original undeformed thickness. Figure 2 depicts the geometry of these quasi-spherical lesions, including r(t) and h(t).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 h:
(1)w=Hλ(t)2(c1log[1-(eE11-1c2)]+c3E11),
where c1 and c3 are constants with dimensions N·m−2, c2 is dimensionless, and E11=(λ(t)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 E11=0, yielding
(2)w=(-c1Hc2+c3H)E11+(-c1H2c22+3c1H2c2-2c3H)E112=+O[E113].
If the 2D deformation gradient F is diag(λ(t),λ(t)), then the Cauchy stress resultant tensor T for inner membrane stress can be represented by the following equation [11]:
(3)Tαβ=1detFFαγFβδ∂w∂Eγδ+2μmHλ(t)3dλ(t)dt,α,β,γ,δ=1,2.
Thus, the isotropic stress resultant T (T11) is given by the following:
(4)T=-c1Hc2+c3H+(-c1H2c22+3c1H2c2-2c3H)(λ(t)2-1)=+2μmHλ(t)3dλ(t)dt.


#### 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 Q is
(5)Q=∑i=17Aicos(2πfit+ϕi).
Table 1 contains the constants Ai, fi, and ϕi [12].Table 1
Lung Fourier series parameter table.
i
A (L·s−1)
f (Hz)
ϕ
1
1
0.156250
4.95
2
0.3
0.390625
3.82
3
0.25
0.859375
4.37
4
0.25
1.484375
3.67
5
0.175
2.421875
4.05
6
0.175
4.609375
4.13
7
0.125
8.046875
4.02Note 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ΔP=RbQ, where Rb was the respiratory resistance of the bronchiole network. Because Rb is 1.0 mmHg·s·L−1 [13] and the original pressure was atmospheric pressure (760 mmHg), the final pressure series was as follows:
(6)Pinner=760+∑i=17Aicos(2πfit+ϕi).


#### 2.1.3. Dynamic Equation System

We obtained the dynamic equation for the wall through linear momentum balance to yield(7)ρmhRd2λ(t)dt2=P(t)-2T(λ(t))λ(t)R,
where ρm is the density of the membrane and P(t) 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]:
(8)Pouter=p∞+ρfR2(λ(t)d2λ(t)dt2+32(dλ(t)dt)2),
where p∞ is a constant value representing the pressure at a significantly large distance from the lung and ρf 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 μf, and viscoelasticity, we arrived at the following differential equation for the air-tissue-pleural fluid system:
(9)(ρmHRλ(t)2+ρfR2λ(t))d2λ(t)dt2+32ρfR2(dλ(t)dt)2+4μfλ(t)×dλ(t)dt+2T(λ(t))λ(t)R+4μmHλ(t)4Rdλ(t)dt=Pinner(t)-p∞,
where μf and μm 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:
(10)y0≡λ(t),τ≡t|c1|HρmR2H,b≡ρfRρmH,d≡HR,mf≡μfρm|c1|H2,m≡μmρm|c1|H2,F(τ)≡(Pinner(t)-p∞)Rc1H,f≡Tc1H.
The resulting nondimensionalized equation system was as follows:
(11)dy0dτ=y1,dy1dτ=F(τ)-3by12/2-4mfy1/y0-2f(y0)/y0-4mdy1/y04by0+y0-2.


#### 2.1.5. Parameter Table

Table2 is a complete collection of the relevant parameter values for the model.Table 2
List of parameters.
Parameter (variable)
Value
Citation
Material parameters (c1, c2, c3)
-
22.5
×
1
0
5 N·m−2, 1.26, -7.8×105 N·m−2
Denny and Schroter, 2006 [8]
Undeformed lesion thickness (H)
10
-
3 m
Amjadi et al., 2007 [15]
Undeformed lesion radius (R)
10
-
2 m
Amjadi et al., 2007 [15]
Resistance of bronchiole network (Rb)
1 mmHg·s·L−1
Ben-Tal, 2006 [13]
Pressure at infinity (p∞)
0 mmHg
Imposed
Pleural fluid viscosity (μf)
1.39–1.57×10-3 Pa·s
Yetkin et al., 2007 [16]
Pleural fluid density (ρf)
980.3 kg·m−3
Rubins and Manning, “Pleural Effusion Workup” [17]
Solid membrane viscosity (μm)
7
×
1
0
-
2 Pa·s
Girnyk et al., 2006 [18]
Membrane density (ρm)
1050 kg·m−2
Shah and Humphrey, 1999 [14]
Atmospheric pressure
760 mmHg
Ben-Tal, 2006 [13]


### 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λ(t) versus time reveal sustained stable oscillations of the membrane for normal lung parenchymal tissue parameters of the Denny-Schroter model: c1=-22.5×105, c2=1.26, and c3=-7.8×105. 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 c2=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 c2∈[.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 c2≈0.78 and then decreases monotonically from that value.Figure 3
Stretch factorλ(t) for different stiffness constant (c2). This figure shows three plots of stretch factor with respect to dimensionless time for the stiffness constants c2=1.00,1.26,and1.52. Variations in the stiffness constant clearly affect the amplitude and frequency of oscillations.Quantified effects of stiffness constant (c2) 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 c2≈0.78.
(a)
(b)Figure5 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 λ(t) versus λ˙(t) further suggests the existence of an orbit about a stable equilibria. The spiral nature of the superimposed vector field, representing the stretch rate λ˙(t) versus the derivative of the stretch rate λ¨(t), 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 c1=-22.5×105, c2=1.26, and c3=-1.39×106, but a rigorous mathematical analysis had been called for to fully determine whether the instabilities were truly physical or just due to numerical effects.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 Section2.4 for more information about physically unstable situations.
(a)
(b)Relative force proportion graph. (a) The graph depicts the relative magnitude of the dimensionless forces at play over dimensionless time: pressure(F(τ)), radial stress (4mfy1/y0), internal membrane (2f(y0)/y0), fluid structure (3by12/2), and viscoelastic (4mdy1/y0). 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.
(a)
(b)


### 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 approximatingF(t)≈F=RP/c1H. The system is at an equilibrium for any point (λ(t),λ˙(t)) such that y˙0=0 and y˙1=0. This implies that λ˙(t)=0 for any equilibrium point. The equilibrium stretch ratio λeq may be found by solving F-(f(λ(t))/λ(t))=0 as follows:
(12)λeq=-α±-4βγ+α22γ,
where α=c1c22F, β=-c1+5c1c2-6c22c3, and γ=c1-3c1c2+4c22c3. For physically meaningful values of λ(t), 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 matrixJ, ODE instability is prescribed for detJ<0ortrJ>0. For the previous dynamical system in y0 and y1, the Jacobian matrix about the positive equilibrium point reduces to(13)J=[012Fγ(-4βγ+α2-α-4βγ+α2)α(-α+-4βγ+α2)2((4γ2/(-α+-4βγ+α2)2)+(b(-α+-4βγ+α2)/2γ))(2mf(α+-4βγ+α2)/β)-(64γ4dm/(-α+-4βγ+α2)4)((4γ2/(-α+-4βγ+α2)2)+(b(-α+-4βγ+α2)/2γ))].The constraint of physically relevant material parameters imposes the inequalitiesc1<0,  c2>0, and c3<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 F=AP/c1H, α is always greater than zero for physically relevant material constants (A,H,P>0).Theorem 1.
For any values ofγ and β such that γ<0,  β>0, and γ<b1/3α-b2/3β, the system is unstable because detJ<0.The expressiondetJ=-∂y˙1/∂y0 is less than zero for all β,γ such that sgn(∂y˙1/∂y0)>0. The sign of this expression is in turn governed by the following three subexpressions. The subexpression of (-α+-4βγ+α2)2 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 -c1+5c1c2-6c22c3>0. Isolation of c3 reduces the inequality to c3<c1(5c2-1)/6c22. Likewise, γ=c1-3c1c2+4c22c3>0 reduces to c3>c1(3c2-1)/4c22 and places bounds on c3. However, the inequality c1(3c2-1)/4c22<c1(5c2-1)/6c22 reduces to 0<c1(c2+1), which cannot be true from the restrictions of c1<0 and c2>0.Lemma 3.
The expression(14)-4βγ+α2-α-4βγ+α2
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βγ+α2-αα2-4βγ>0 reduces to α2-4βγ>αα2-4βγ. Because both sides of the inequality are positive, the inequality may be squared and the terms brought to one side to give (α2-4βγ)(-4βγ)>0. Because only values of α,β,   and γ such that α2-4βγ>0 are considered, the expression is true βγ<0. Thus, -4βγ+α2-α-4βγ+α2>0 for either β<0 and γ>0 or β>0 and γ<0. Likewise, -4βγ+α2-α-4βγ+α2<0 for β<0 and γ<0 or β>0 and γ>0.Lemma 4.
For physically relevant restriction ofb>0, the expression
(15)4γ2(-α+-4βγ+α2)2+b(-α+-4βγ+α2)2γ
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; γ<b1/3α-b2/3β,  γ<0, and β>0.
For either the conditions ofγ>0 and β<0 or γ<0 and β<0, the expression b(-α+-4βγ+α2)/2γ is always greater than zero. Hence, because 4γ2/(-α+-4βγ+α2)2 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γ3+b(-α+-4βγ+α2)<0 or γ<-(b1/3/2)(-α+-4βγ+α2). This in turn reduces to 4γ/b2/3-4α/b1/3+4β<0 or γ<b1/3α-b2/3β.Lemma 5.
For physically relevant restriction ofb>0, the expression (15) is less than zero for any value of β,γ such that γ<0,  β>0, and γ>b1/3α-b2/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 γ>b1/3α-b2/3β is proven likewise as in Lemma 4.Lemma 6.
The expression(16)2Fγα
is greater than zero for all γ<0 and less than zero for all γ>0. By definition, F<0 and α>0, and thus the result follows.
The expressiondetJ<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 γ<b1/3α-b2/3β the system is unstable because detJ<0.Theorem 7.
For any values ofγ and β such that γ<0, β>0, γ>b1/3α-b2/3β, and β>αγ2(dm/mf)1/3+γ(dm/mf)2/3, the system is unstable because trJ>0.The expressiontrJ=-∂y˙1/∂y1 is less than zero for all β,γ such that sgn(∂y˙1/∂y1)>0. The sign of this expression is in turn governed by subexpressions (15) and (17).Lemma 8.
For physically relevant parameters (mf,m,d>0),
(17)2mf(α+-4βγ+α2)β-64γ4dm(-α+-4βγ+α2)4
is less than zero for all β<0 or for β>0,  γ<0, and β>αγ2(dm/mf)1/3+γ(dm/mf)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
(18)2mfβ(α+α2-4βγ)(-α+α2-4βγ)4<64γ4dm
or (-α+α2-4βγ)3>-8γ3dm/mf. This in turn yields α2-4βγ>(α-2γ(dm/mf)1/3)2, which upon isolation of β gives the expression β>αγ2(dm/mf)1/3+γ(dm/mf)2/3.The expressiontrJ>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, γ>b1/3α-b2/3β, and β>αγ2(dm/mf)1/3+γ(dm/mf)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 c1, c2, and c3.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].
(a)
(b)
(c)


### 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Φ(Jm-3). Similar to Volokh and Vorp [19], we built a finite energy constitutive model Ψ based on the framework of the hyperelastic model w which was constructed as follows:
(19)Ψ=HΦ(1-exp(-wHΦ)).
Note that as w→∞, Ψ→HΦ a finite value of energy before rupture in contrast to the hyperelastic model where w may grow unbounded as λ(t)→∞. Similarly, in the limiting case HΦ≪w,  Ψ≈w as can be seen via Taylor expansion. Reformulating the Cauchy stress resultant tensor T given by [14] with the new strain-energy function yields
(20)Tαβ′=1detFFαγFβδ∂Ψ∂Eγδ+2μmHλ(t)3dλ(t)dt,α,β,γ,δ=1,2.
Thus, the new isotropic stress resultant T′ (T11′) is given by the following:
(21)T′=(-c1Hc2+c3H+(-c1H2c22+3c1H2c2-2c3H)(λ(t)2-1))=×exp((c1H/2c22-3c1H/2c2+2c3H)((λ(t)2-1)2/4)HΦ(c1H/c2-c3H)((λ(t)2-1)/2)HΦ×expmM+(c1H/2c22-3c1H/2c2+2c3H)((λ(t)2-1)2/4)HΦ)=+2μmHλ(t)3dλ(t)dt.As demonstrated in Figure8, 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(c1,c2,c3) 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:
(22)w′=Ac(τ)(Hλ(t)2(c1log[1-(eE11-1c2)]+c3E11))w′=+Ae(τ)(bHλ(t)2(λ(t)-logλ(t)-1)),w′≈(-Ac(τ)c1Hc2+Ac(τ)c3H)λ(t)2-12w′≈+(Ae(τ)bH2-Ac(τ)c1H2c22+3Ac(τ)c1H2c2w′≈+M-2Ac(τ)c3HAe(τ)bH2)(λ(t)2-1)24,
where Ac(τ) and Ae(τ) 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 C and elastin E defined by ρcAc and ρeAe (where ρc and ρe are equal to 5.10 kg·m−3 [21]), the new momentum balance was as follows:
(23)ρmhRd2λ(t)dt2+d(Aeρe+Acρc)dtvrH=P(t)-2T(λ(t))λ(t)R,
where vr(t) is the velocity at the boundary of an expanding biological membrane given by Rdλ(t)/dt. With the nondimensionalized term mp=(A/ρm|c1|)(dC/dt+dE/dt), the first two differential equations of the dynamical system became
(24)dy0dτ=y1,dy1dτ=F(τ)-3by02/2-4mfy1/y0-2f(y0)/y0-4mdy1/y04-mpy1/y02by0+y0-2.
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)dAcdt=k1Ac|ϕ′-ψ′|n,dAedt=-k2Ae|ϕ′-ψ′|n.
Because k1 and k2 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.Collagen-elastin dynamics may lead to escapes from stable orbits. (a) The graph showsλ(τ) with the additional collagen-elastin dynamics calibrated with Ac=100,  Ae=50,  k1=1×107,  k2=1×103,  and n=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.
(a)
(b)


### 2.6. Parameter Estimation fromIn Vivo Data

In order to build a database for parameter estimation on the proposed model, we collectedin 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, Gx and Gy, through application of a matrix mask in order to detect edges. Pixels that meet a certain threshold for G  (|Gx|+|Gy|) 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)ρmv2H), 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 parametersc1, c2, and c3 as well as the corresponding coefficients of determination (R2 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.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.
(a)
(b)
(c)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.


## 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λ(t) was defined as r(t)/R, where r(t) represents the deformed radius varying with time and R is the original undeformed radius. Assuming membrane incompressibility, the constant volume condition V=4πr2(t)h(t) implies that the deformed thickness is h(t)=H/λ(t)2, where H represents the original undeformed thickness. Figure 2 depicts the geometry of these quasi-spherical lesions, including r(t) and h(t).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 h:
(1)w=Hλ(t)2(c1log[1-(eE11-1c2)]+c3E11),
where c1 and c3 are constants with dimensions N·m−2, c2 is dimensionless, and E11=(λ(t)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 E11=0, yielding
(2)w=(-c1Hc2+c3H)E11+(-c1H2c22+3c1H2c2-2c3H)E112=+O[E113].
If the 2D deformation gradient F is diag(λ(t),λ(t)), then the Cauchy stress resultant tensor T for inner membrane stress can be represented by the following equation [11]:
(3)Tαβ=1detFFαγFβδ∂w∂Eγδ+2μmHλ(t)3dλ(t)dt,α,β,γ,δ=1,2.
Thus, the isotropic stress resultant T (T11) is given by the following:
(4)T=-c1Hc2+c3H+(-c1H2c22+3c1H2c2-2c3H)(λ(t)2-1)=+2μmHλ(t)3dλ(t)dt.


### 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 Q is
(5)Q=∑i=17Aicos(2πfit+ϕi).
Table 1 contains the constants Ai, fi, and ϕi [12].Table 1
Lung Fourier series parameter table.
i
A (L·s−1)
f (Hz)
ϕ
1
1
0.156250
4.95
2
0.3
0.390625
3.82
3
0.25
0.859375
4.37
4
0.25
1.484375
3.67
5
0.175
2.421875
4.05
6
0.175
4.609375
4.13
7
0.125
8.046875
4.02Note 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ΔP=RbQ, where Rb was the respiratory resistance of the bronchiole network. Because Rb is 1.0 mmHg·s·L−1 [13] and the original pressure was atmospheric pressure (760 mmHg), the final pressure series was as follows:
(6)Pinner=760+∑i=17Aicos(2πfit+ϕi).


### 2.1.3. Dynamic Equation System

We obtained the dynamic equation for the wall through linear momentum balance to yield(7)ρmhRd2λ(t)dt2=P(t)-2T(λ(t))λ(t)R,
where ρm is the density of the membrane and P(t) 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]:
(8)Pouter=p∞+ρfR2(λ(t)d2λ(t)dt2+32(dλ(t)dt)2),
where p∞ is a constant value representing the pressure at a significantly large distance from the lung and ρf 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 μf, and viscoelasticity, we arrived at the following differential equation for the air-tissue-pleural fluid system:
(9)(ρmHRλ(t)2+ρfR2λ(t))d2λ(t)dt2+32ρfR2(dλ(t)dt)2+4μfλ(t)×dλ(t)dt+2T(λ(t))λ(t)R+4μmHλ(t)4Rdλ(t)dt=Pinner(t)-p∞,
where μf and μm 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:
(10)y0≡λ(t),τ≡t|c1|HρmR2H,b≡ρfRρmH,d≡HR,mf≡μfρm|c1|H2,m≡μmρm|c1|H2,F(τ)≡(Pinner(t)-p∞)Rc1H,f≡Tc1H.
The resulting nondimensionalized equation system was as follows:
(11)dy0dτ=y1,dy1dτ=F(τ)-3by12/2-4mfy1/y0-2f(y0)/y0-4mdy1/y04by0+y0-2.


### 2.1.5. Parameter Table

Table2 is a complete collection of the relevant parameter values for the model.Table 2
List of parameters.
Parameter (variable)
Value
Citation
Material parameters (c1, c2, c3)
-
22.5
×
1
0
5 N·m−2, 1.26, -7.8×105 N·m−2
Denny and Schroter, 2006 [8]
Undeformed lesion thickness (H)
10
-
3 m
Amjadi et al., 2007 [15]
Undeformed lesion radius (R)
10
-
2 m
Amjadi et al., 2007 [15]
Resistance of bronchiole network (Rb)
1 mmHg·s·L−1
Ben-Tal, 2006 [13]
Pressure at infinity (p∞)
0 mmHg
Imposed
Pleural fluid viscosity (μf)
1.39–1.57×10-3 Pa·s
Yetkin et al., 2007 [16]
Pleural fluid density (ρf)
980.3 kg·m−3
Rubins and Manning, “Pleural Effusion Workup” [17]
Solid membrane viscosity (μm)
7
×
1
0
-
2 Pa·s
Girnyk et al., 2006 [18]
Membrane density (ρm)
1050 kg·m−2
Shah and Humphrey, 1999 [14]
Atmospheric pressure
760 mmHg
Ben-Tal, 2006 [13]


## 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 h:
(1)w=Hλ(t)2(c1log[1-(eE11-1c2)]+c3E11),
where c1 and c3 are constants with dimensions N·m−2, c2 is dimensionless, and E11=(λ(t)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 E11=0, yielding
(2)w=(-c1Hc2+c3H)E11+(-c1H2c22+3c1H2c2-2c3H)E112=+O[E113].
If the 2D deformation gradient F is diag(λ(t),λ(t)), then the Cauchy stress resultant tensor T for inner membrane stress can be represented by the following equation [11]:
(3)Tαβ=1detFFαγFβδ∂w∂Eγδ+2μmHλ(t)3dλ(t)dt,α,β,γ,δ=1,2.
Thus, the isotropic stress resultant T (T11) is given by the following:
(4)T=-c1Hc2+c3H+(-c1H2c22+3c1H2c2-2c3H)(λ(t)2-1)=+2μmHλ(t)3dλ(t)dt.


## 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 Q is
(5)Q=∑i=17Aicos(2πfit+ϕi).
Table 1 contains the constants Ai, fi, and ϕi [12].Table 1
Lung Fourier series parameter table.
i
A (L·s−1)
f (Hz)
ϕ
1
1
0.156250
4.95
2
0.3
0.390625
3.82
3
0.25
0.859375
4.37
4
0.25
1.484375
3.67
5
0.175
2.421875
4.05
6
0.175
4.609375
4.13
7
0.125
8.046875
4.02Note 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ΔP=RbQ, where Rb was the respiratory resistance of the bronchiole network. Because Rb is 1.0 mmHg·s·L−1 [13] and the original pressure was atmospheric pressure (760 mmHg), the final pressure series was as follows:
(6)Pinner=760+∑i=17Aicos(2πfit+ϕi).


## 2.1.3. Dynamic Equation System

We obtained the dynamic equation for the wall through linear momentum balance to yield(7)ρmhRd2λ(t)dt2=P(t)-2T(λ(t))λ(t)R,
where ρm is the density of the membrane and P(t) 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]:
(8)Pouter=p∞+ρfR2(λ(t)d2λ(t)dt2+32(dλ(t)dt)2),
where p∞ is a constant value representing the pressure at a significantly large distance from the lung and ρf 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 μf, and viscoelasticity, we arrived at the following differential equation for the air-tissue-pleural fluid system:
(9)(ρmHRλ(t)2+ρfR2λ(t))d2λ(t)dt2+32ρfR2(dλ(t)dt)2+4μfλ(t)×dλ(t)dt+2T(λ(t))λ(t)R+4μmHλ(t)4Rdλ(t)dt=Pinner(t)-p∞,
where μf and μm 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:
(10)y0≡λ(t),τ≡t|c1|HρmR2H,b≡ρfRρmH,d≡HR,mf≡μfρm|c1|H2,m≡μmρm|c1|H2,F(τ)≡(Pinner(t)-p∞)Rc1H,f≡Tc1H.
The resulting nondimensionalized equation system was as follows:
(11)dy0dτ=y1,dy1dτ=F(τ)-3by12/2-4mfy1/y0-2f(y0)/y0-4mdy1/y04by0+y0-2.


## 2.1.5. Parameter Table

Table2 is a complete collection of the relevant parameter values for the model.Table 2
List of parameters.
Parameter (variable)
Value
Citation
Material parameters (c1, c2, c3)
-
22.5
×
1
0
5 N·m−2, 1.26, -7.8×105 N·m−2
Denny and Schroter, 2006 [8]
Undeformed lesion thickness (H)
10
-
3 m
Amjadi et al., 2007 [15]
Undeformed lesion radius (R)
10
-
2 m
Amjadi et al., 2007 [15]
Resistance of bronchiole network (Rb)
1 mmHg·s·L−1
Ben-Tal, 2006 [13]
Pressure at infinity (p∞)
0 mmHg
Imposed
Pleural fluid viscosity (μf)
1.39–1.57×10-3 Pa·s
Yetkin et al., 2007 [16]
Pleural fluid density (ρf)
980.3 kg·m−3
Rubins and Manning, “Pleural Effusion Workup” [17]
Solid membrane viscosity (μm)
7
×
1
0
-
2 Pa·s
Girnyk et al., 2006 [18]
Membrane density (ρm)
1050 kg·m−2
Shah and Humphrey, 1999 [14]
Atmospheric pressure
760 mmHg
Ben-Tal, 2006 [13]


## 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λ(t) versus time reveal sustained stable oscillations of the membrane for normal lung parenchymal tissue parameters of the Denny-Schroter model: c1=-22.5×105, c2=1.26, and c3=-7.8×105. 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 c2=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 c2∈[.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 c2≈0.78 and then decreases monotonically from that value.Figure 3
Stretch factorλ(t) for different stiffness constant (c2). This figure shows three plots of stretch factor with respect to dimensionless time for the stiffness constants c2=1.00,1.26,and1.52. Variations in the stiffness constant clearly affect the amplitude and frequency of oscillations.Quantified effects of stiffness constant (c2) 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 c2≈0.78.
(a)
(b)Figure5 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 λ(t) versus λ˙(t) further suggests the existence of an orbit about a stable equilibria. The spiral nature of the superimposed vector field, representing the stretch rate λ˙(t) versus the derivative of the stretch rate λ¨(t), 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 c1=-22.5×105, c2=1.26, and c3=-1.39×106, but a rigorous mathematical analysis had been called for to fully determine whether the instabilities were truly physical or just due to numerical effects.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 Section2.4 for more information about physically unstable situations.
(a)
(b)Relative force proportion graph. (a) The graph depicts the relative magnitude of the dimensionless forces at play over dimensionless time: pressure(F(τ)), radial stress (4mfy1/y0), internal membrane (2f(y0)/y0), fluid structure (3by12/2), and viscoelastic (4mdy1/y0). 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.
(a)
(b)


## 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 approximatingF(t)≈F=RP/c1H. The system is at an equilibrium for any point (λ(t),λ˙(t)) such that y˙0=0 and y˙1=0. This implies that λ˙(t)=0 for any equilibrium point. The equilibrium stretch ratio λeq may be found by solving F-(f(λ(t))/λ(t))=0 as follows:
(12)λeq=-α±-4βγ+α22γ,
where α=c1c22F, β=-c1+5c1c2-6c22c3, and γ=c1-3c1c2+4c22c3. For physically meaningful values of λ(t), 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 matrixJ, ODE instability is prescribed for detJ<0ortrJ>0. For the previous dynamical system in y0 and y1, the Jacobian matrix about the positive equilibrium point reduces to(13)J=[012Fγ(-4βγ+α2-α-4βγ+α2)α(-α+-4βγ+α2)2((4γ2/(-α+-4βγ+α2)2)+(b(-α+-4βγ+α2)/2γ))(2mf(α+-4βγ+α2)/β)-(64γ4dm/(-α+-4βγ+α2)4)((4γ2/(-α+-4βγ+α2)2)+(b(-α+-4βγ+α2)/2γ))].The constraint of physically relevant material parameters imposes the inequalitiesc1<0,  c2>0, and c3<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 F=AP/c1H, α is always greater than zero for physically relevant material constants (A,H,P>0).Theorem 1.
For any values ofγ and β such that γ<0,  β>0, and γ<b1/3α-b2/3β, the system is unstable because detJ<0.The expressiondetJ=-∂y˙1/∂y0 is less than zero for all β,γ such that sgn(∂y˙1/∂y0)>0. The sign of this expression is in turn governed by the following three subexpressions. The subexpression of (-α+-4βγ+α2)2 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 -c1+5c1c2-6c22c3>0. Isolation of c3 reduces the inequality to c3<c1(5c2-1)/6c22. Likewise, γ=c1-3c1c2+4c22c3>0 reduces to c3>c1(3c2-1)/4c22 and places bounds on c3. However, the inequality c1(3c2-1)/4c22<c1(5c2-1)/6c22 reduces to 0<c1(c2+1), which cannot be true from the restrictions of c1<0 and c2>0.Lemma 3.
The expression(14)-4βγ+α2-α-4βγ+α2
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βγ+α2-αα2-4βγ>0 reduces to α2-4βγ>αα2-4βγ. Because both sides of the inequality are positive, the inequality may be squared and the terms brought to one side to give (α2-4βγ)(-4βγ)>0. Because only values of α,β,   and γ such that α2-4βγ>0 are considered, the expression is true βγ<0. Thus, -4βγ+α2-α-4βγ+α2>0 for either β<0 and γ>0 or β>0 and γ<0. Likewise, -4βγ+α2-α-4βγ+α2<0 for β<0 and γ<0 or β>0 and γ>0.Lemma 4.
For physically relevant restriction ofb>0, the expression
(15)4γ2(-α+-4βγ+α2)2+b(-α+-4βγ+α2)2γ
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; γ<b1/3α-b2/3β,  γ<0, and β>0.
For either the conditions ofγ>0 and β<0 or γ<0 and β<0, the expression b(-α+-4βγ+α2)/2γ is always greater than zero. Hence, because 4γ2/(-α+-4βγ+α2)2 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γ3+b(-α+-4βγ+α2)<0 or γ<-(b1/3/2)(-α+-4βγ+α2). This in turn reduces to 4γ/b2/3-4α/b1/3+4β<0 or γ<b1/3α-b2/3β.Lemma 5.
For physically relevant restriction ofb>0, the expression (15) is less than zero for any value of β,γ such that γ<0,  β>0, and γ>b1/3α-b2/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 γ>b1/3α-b2/3β is proven likewise as in Lemma 4.Lemma 6.
The expression(16)2Fγα
is greater than zero for all γ<0 and less than zero for all γ>0. By definition, F<0 and α>0, and thus the result follows.
The expressiondetJ<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 γ<b1/3α-b2/3β the system is unstable because detJ<0.Theorem 7.
For any values ofγ and β such that γ<0, β>0, γ>b1/3α-b2/3β, and β>αγ2(dm/mf)1/3+γ(dm/mf)2/3, the system is unstable because trJ>0.The expressiontrJ=-∂y˙1/∂y1 is less than zero for all β,γ such that sgn(∂y˙1/∂y1)>0. The sign of this expression is in turn governed by subexpressions (15) and (17).Lemma 8.
For physically relevant parameters (mf,m,d>0),
(17)2mf(α+-4βγ+α2)β-64γ4dm(-α+-4βγ+α2)4
is less than zero for all β<0 or for β>0,  γ<0, and β>αγ2(dm/mf)1/3+γ(dm/mf)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
(18)2mfβ(α+α2-4βγ)(-α+α2-4βγ)4<64γ4dm
or (-α+α2-4βγ)3>-8γ3dm/mf. This in turn yields α2-4βγ>(α-2γ(dm/mf)1/3)2, which upon isolation of β gives the expression β>αγ2(dm/mf)1/3+γ(dm/mf)2/3.The expressiontrJ>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, γ>b1/3α-b2/3β, and β>αγ2(dm/mf)1/3+γ(dm/mf)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 c1, c2, and c3.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].
(a)
(b)
(c)


## 2.3.1. Determining System Equilibria

In order to simplify the following equilibria identification and stability analysis, we made the system autonomous by approximatingF(t)≈F=RP/c1H. The system is at an equilibrium for any point (λ(t),λ˙(t)) such that y˙0=0 and y˙1=0. This implies that λ˙(t)=0 for any equilibrium point. The equilibrium stretch ratio λeq may be found by solving F-(f(λ(t))/λ(t))=0 as follows:
(12)λeq=-α±-4βγ+α22γ,
where α=c1c22F, β=-c1+5c1c2-6c22c3, and γ=c1-3c1c2+4c22c3. For physically meaningful values of λ(t), 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 matrixJ, ODE instability is prescribed for detJ<0ortrJ>0. For the previous dynamical system in y0 and y1, the Jacobian matrix about the positive equilibrium point reduces to(13)J=[012Fγ(-4βγ+α2-α-4βγ+α2)α(-α+-4βγ+α2)2((4γ2/(-α+-4βγ+α2)2)+(b(-α+-4βγ+α2)/2γ))(2mf(α+-4βγ+α2)/β)-(64γ4dm/(-α+-4βγ+α2)4)((4γ2/(-α+-4βγ+α2)2)+(b(-α+-4βγ+α2)/2γ))].The constraint of physically relevant material parameters imposes the inequalitiesc1<0,  c2>0, and c3<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 F=AP/c1H, α is always greater than zero for physically relevant material constants (A,H,P>0).Theorem 1.
For any values ofγ and β such that γ<0,  β>0, and γ<b1/3α-b2/3β, the system is unstable because detJ<0.The expressiondetJ=-∂y˙1/∂y0 is less than zero for all β,γ such that sgn(∂y˙1/∂y0)>0. The sign of this expression is in turn governed by the following three subexpressions. The subexpression of (-α+-4βγ+α2)2 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 -c1+5c1c2-6c22c3>0. Isolation of c3 reduces the inequality to c3<c1(5c2-1)/6c22. Likewise, γ=c1-3c1c2+4c22c3>0 reduces to c3>c1(3c2-1)/4c22 and places bounds on c3. However, the inequality c1(3c2-1)/4c22<c1(5c2-1)/6c22 reduces to 0<c1(c2+1), which cannot be true from the restrictions of c1<0 and c2>0.Lemma 3.
The expression(14)-4βγ+α2-α-4βγ+α2
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βγ+α2-αα2-4βγ>0 reduces to α2-4βγ>αα2-4βγ. Because both sides of the inequality are positive, the inequality may be squared and the terms brought to one side to give (α2-4βγ)(-4βγ)>0. Because only values of α,β,   and γ such that α2-4βγ>0 are considered, the expression is true βγ<0. Thus, -4βγ+α2-α-4βγ+α2>0 for either β<0 and γ>0 or β>0 and γ<0. Likewise, -4βγ+α2-α-4βγ+α2<0 for β<0 and γ<0 or β>0 and γ>0.Lemma 4.
For physically relevant restriction ofb>0, the expression
(15)4γ2(-α+-4βγ+α2)2+b(-α+-4βγ+α2)2γ
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; γ<b1/3α-b2/3β,  γ<0, and β>0.
For either the conditions ofγ>0 and β<0 or γ<0 and β<0, the expression b(-α+-4βγ+α2)/2γ is always greater than zero. Hence, because 4γ2/(-α+-4βγ+α2)2 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γ3+b(-α+-4βγ+α2)<0 or γ<-(b1/3/2)(-α+-4βγ+α2). This in turn reduces to 4γ/b2/3-4α/b1/3+4β<0 or γ<b1/3α-b2/3β.Lemma 5.
For physically relevant restriction ofb>0, the expression (15) is less than zero for any value of β,γ such that γ<0,  β>0, and γ>b1/3α-b2/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 γ>b1/3α-b2/3β is proven likewise as in Lemma 4.Lemma 6.
The expression(16)2Fγα
is greater than zero for all γ<0 and less than zero for all γ>0. By definition, F<0 and α>0, and thus the result follows.
The expressiondetJ<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 γ<b1/3α-b2/3β the system is unstable because detJ<0.Theorem 7.
For any values ofγ and β such that γ<0, β>0, γ>b1/3α-b2/3β, and β>αγ2(dm/mf)1/3+γ(dm/mf)2/3, the system is unstable because trJ>0.The expressiontrJ=-∂y˙1/∂y1 is less than zero for all β,γ such that sgn(∂y˙1/∂y1)>0. The sign of this expression is in turn governed by subexpressions (15) and (17).Lemma 8.
For physically relevant parameters (mf,m,d>0),
(17)2mf(α+-4βγ+α2)β-64γ4dm(-α+-4βγ+α2)4
is less than zero for all β<0 or for β>0,  γ<0, and β>αγ2(dm/mf)1/3+γ(dm/mf)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
(18)2mfβ(α+α2-4βγ)(-α+α2-4βγ)4<64γ4dm
or (-α+α2-4βγ)3>-8γ3dm/mf. This in turn yields α2-4βγ>(α-2γ(dm/mf)1/3)2, which upon isolation of β gives the expression β>αγ2(dm/mf)1/3+γ(dm/mf)2/3.The expressiontrJ>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, γ>b1/3α-b2/3β, and β>αγ2(dm/mf)1/3+γ(dm/mf)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 c1, c2, and c3.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].
(a)
(b)
(c)


## 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Φ(Jm-3). Similar to Volokh and Vorp [19], we built a finite energy constitutive model Ψ based on the framework of the hyperelastic model w which was constructed as follows:
(19)Ψ=HΦ(1-exp(-wHΦ)).
Note that as w→∞, Ψ→HΦ a finite value of energy before rupture in contrast to the hyperelastic model where w may grow unbounded as λ(t)→∞. Similarly, in the limiting case HΦ≪w,  Ψ≈w as can be seen via Taylor expansion. Reformulating the Cauchy stress resultant tensor T given by [14] with the new strain-energy function yields
(20)Tαβ′=1detFFαγFβδ∂Ψ∂Eγδ+2μmHλ(t)3dλ(t)dt,α,β,γ,δ=1,2.
Thus, the new isotropic stress resultant T′ (T11′) is given by the following:
(21)T′=(-c1Hc2+c3H+(-c1H2c22+3c1H2c2-2c3H)(λ(t)2-1))=×exp((c1H/2c22-3c1H/2c2+2c3H)((λ(t)2-1)2/4)HΦ(c1H/c2-c3H)((λ(t)2-1)/2)HΦ×expmM+(c1H/2c22-3c1H/2c2+2c3H)((λ(t)2-1)2/4)HΦ)=+2μmHλ(t)3dλ(t)dt.As demonstrated in Figure8, 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(c1,c2,c3) 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:
(22)w′=Ac(τ)(Hλ(t)2(c1log[1-(eE11-1c2)]+c3E11))w′=+Ae(τ)(bHλ(t)2(λ(t)-logλ(t)-1)),w′≈(-Ac(τ)c1Hc2+Ac(τ)c3H)λ(t)2-12w′≈+(Ae(τ)bH2-Ac(τ)c1H2c22+3Ac(τ)c1H2c2w′≈+M-2Ac(τ)c3HAe(τ)bH2)(λ(t)2-1)24,
where Ac(τ) and Ae(τ) 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 C and elastin E defined by ρcAc and ρeAe (where ρc and ρe are equal to 5.10 kg·m−3 [21]), the new momentum balance was as follows:
(23)ρmhRd2λ(t)dt2+d(Aeρe+Acρc)dtvrH=P(t)-2T(λ(t))λ(t)R,
where vr(t) is the velocity at the boundary of an expanding biological membrane given by Rdλ(t)/dt. With the nondimensionalized term mp=(A/ρm|c1|)(dC/dt+dE/dt), the first two differential equations of the dynamical system became
(24)dy0dτ=y1,dy1dτ=F(τ)-3by02/2-4mfy1/y0-2f(y0)/y0-4mdy1/y04-mpy1/y02by0+y0-2.
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)dAcdt=k1Ac|ϕ′-ψ′|n,dAedt=-k2Ae|ϕ′-ψ′|n.
Because k1 and k2 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.Collagen-elastin dynamics may lead to escapes from stable orbits. (a) The graph showsλ(τ) with the additional collagen-elastin dynamics calibrated with Ac=100,  Ae=50,  k1=1×107,  k2=1×103,  and n=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.
(a)
(b)


## 2.6. Parameter Estimation fromIn Vivo Data

In order to build a database for parameter estimation on the proposed model, we collectedin 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, Gx and Gy, through application of a matrix mask in order to detect edges. Pixels that meet a certain threshold for G  (|Gx|+|Gy|) 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)ρmv2H), 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 parametersc1, c2, and c3 as well as the corresponding coefficients of determination (R2 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.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.
(a)
(b)
(c)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.



---

*Source: 679308-2013-05-14.xml*