Airway Wall Deformation Under Smooth Muscle Contraction

The objective of this thesis is the development and evaluation of a mathematical model that captures the deformation behavior of the human airway under smooth muscle contraction. The problem consists of a multilayered circular cylindrical tube subjected to an active contractile stress generated by the smooth muscle layer. The problem formulation assumes axisymmetric deformation for an incompressible isotropic neo-Hookean material. The model is formulated as fully nonlinear and accommodates large deformation, since these are expected in the human airway and other soft biological materials. The equilibrium equations are obtained, as well as the strain-displacement relationship, boundary conditions, and constitutive equations that describe the contractile stress for neo-Hookean hyperelastic solids. The main focus of the project is to obtain the change in airway caliber after a prescribed contractile stress is introduced in the smooth muscle layer. A MATLAB code was written to numerically solve the resulting nonlinear algebraic equation, which obtains the change in airway caliber. The code allows the user to analyze and change the material parameters that govern the solution, as well as input different dimensions for the multilayer cylinder. In the deformation of the airway, the lung parenchyma plays a very important role in preventing full closure of the airway. Here, the lung parenchyma is treated as an infinitely large continuous solid, in which the airway is embedded. The airway generations selected for further investigation are generations 0, 4, 8, 12 and 16. For the lung parenchyma infinite domain assumption, the percent airway caliber change for the critical contractile stresses of the airway generations were 23.86%, 23.68%, 23.95%, 24.56%, and 25.51% for generations 0,4,8,12, and 16, respectively. The mathematical model does not account for the buckling that is observed in the deformation of the human airway, since the deformation is assumed to be axisymmetric. When the contractile stress in the smooth muscle layer exceeds a certain value, the airway buckles and folds. Therefore, finite element simulations were conducted using ABAQUS software to determine the critical contractile stress at which the airway will buckle. Buckling analyses were performed in order to obtain the buckling modes of the airway, as well as the critical load that will cause the it to buckle, for airway generations 0, 4, 8, 12, and 16. The critical load in which the airway buckles and folds for these generations of the airway are 1.755 kPa, 1.572 kPa, 1.493 kPa, 1.426 kPa, and 1.374 kPa respectively. Contractile stresses used in the analytical model that are larger than the critical buckling loads will provide invalid results since the assumptions made to develop the model will no longer be valid. After the buckling analysis, the first buckling mode is introduced in a post buckling analysis in order to obtain the average change in caliber and the active contractile stresses.

. Asthma is a disease of the airway, classified as the obstruction of air flow due to airway thickening from scarring and inflammation. It is also classified as bronchoconstriction, which is the narrowing of the airways in the lungs due to tightening of the surrounding smooth muscle.   [4], which aimed to obtain the elastic properties of the airway wall for the use in theoretical models. In the study by Codd, strips of trachea from sheep were tested in tension, in both the axial and circumferential directions. These studies did not aim to develop a theoretical model of their own. The studies aimed at obtaining elastic properties in order to be used in future theoretical models.
A theoretical model developed by Lai-Fook and his colleagues [5] tried to predict the diameter change of the airway wall. The model separates the airway into two portions, one which examines the elastic behavior of the parenchyma, and the other of the airway wall, separately. Once both have been studied, then the two models can be used to describe the entirety of the system. Initially, the derivation considers a homogeneous parenchymal structure, separated from the airway wall, and treated as an elastic continuum. As the volume within the lung varies, so does the imaginary cylindrical surface that runs through the lung parenchyma. According to Lai-Fook, the diameter of the imaginary hole varies as the cube root of the lung volume. Since lung volume is a function of trans-pulmonary pressure, then it is possible to construct a pressure-diameter relationship for the imaginary hole. When the imaginary hole is replaced by the airway wall, then the change in diameter will also depend on the elastic characteristic of the airway wall.
Other theoretical models were developed from the need to have a tube law for all of the pulmonary airways. Lambert et al. [31] developed a theoretical model based on the pressure-area measurements given by a previous study conducted by Hyatt et al. [7]. Lambert [11] continued the earlier work of the tube law, and fully developed a model which includes a balance between the applied forces from smooth muscle constriction and the structural properties of the lung parenchymal and the airway wall. However, the model heavily relies on the tube law based on sparse extrapolated data.
The theoretical models previously developed mainly deal with pressurevolume relations, or with a tube law derived from extrapolated data. These models do not completely capture the entire deformation and contractile stress of the airway wall.
These models also do not include the complete deformation, buckling, and folding of the airway wall.
Similar to the cylindrical geometry of the human airway, more rigorous theoretical models using continuum mechanics and hyperelastic constitutive relations have been developed to describe this type of problem. Zhu et al [12]  In order to understand the buckling and folding of the airway, analyses were completed using finite element analysis. Lambert [6] performed finite element analyses of the process of buckling with a tube law corresponding to the deformed buckled state. In the model, the inner layer corresponding to the subepithelial collagen layer is modeled as a stiff band. The smooth muscle layer is modeled as a compliant elastic thick layer. The smooth muscle layer is made to constrict uniformly, since the folding of the airway wall depends on the thickness of the stiff and smooth muscle layer, as well as the means that the airway wall constricts. Lambert discovered that for finite element simulations, simple application of external pressure only produced a "peanut shape" buckling, and not the multi-fold buckling that is seen in airways. The multi-lobe buckling is only seen when a band applies a circumferential displacement on the periphery of the smooth muscle layer.
In this study, it is necessary to perform finite element analysis to obtain the critical contractile stress that will cause the airway to buckle. The theoretical model developed will no longer be valid due to the axisymmetric deformation assumption.
Once the airway buckles, the deformation will not be axisymmetric, and as the assumption that the airway caliber is perfectly circular will no longer hold as well.

SECTION 2 -MODEL DESCRIPTION
The physiology of the airway wall is extremely complex. The airway has different sections, as well as layers of tissue that play various roles in breathing. The  A detailed schematic of the human airway is seen in Figure 2   For the rest of this work, the passive load bearing inner-most layer will be referred to as the stiff layer; the contractile middle layer will be referred to as the smooth muscle layer, and the tethering outer layer will be the lung parenchyma layer.
The stiff layer consists of material that is passive in the sense that it will deform and develop internal stress, strain, and displacement in response to internal or external loadings. The smooth muscle layer behaves in a different manner. It is the aim of this project to develop a model that captures the active contractile stress developed by the smooth muscle. It is necessary to separate the smooth muscle layer into two components. One component of the smooth muscle layer is passive and load bearing, similar to that of the stiff layer. The other component applies an active contractile stress that causes the airway to deform and contract. The lung parenchyma layer is also a passive layer, similar to the stiff layer. However, the lung parenchyma layer plays a different role in the deformation of the airway. The constriction of the smooth muscle is due to the force that the muscle can develop in response to a stimulus. In intrapulmonary airway, the decrease in airway caliber due to smooth muscle constriction is opposed by the surrounding lung parenchyma and the stiff layer. This tethering effect of the lung parenchyma on the airway is an elastic constraint that opposes muscle contraction. In the airway, the elastic constraint from the lung parenchyma increases with the increase in lung volume. The shear modulus of the lung parenchyma increases in magnitude with increase in transpulmonary pressure, as well as increase in lung volume [13]. Maximal airway narrowing during constriction is greater at low lung volumes, when the elastic load from the lung parenchyma is lower, than at higher volumes. It can be seen that the material property of the lung parenchyma is a non-linear function of the pressure and lung volume.
It is assumed in this case, when the smooth muscle contract, the deformation of the airway is analyzed under a fixed volume of the lung in order to simplify the nonlinear dependence for the mechanical properties of the lung parenchyma. To obtain the best solution for the effect of the lung parenchyma on airway deformation, two distinct descriptions of the lung parenchyma will be used. In the first formulation of the model, the lung parenchyma will be modeled as an infinite elastic continuum, as seen in figure 2.4. This assumption hinges on the fact that the numbers of parenchymal attachments that surround the airway are large in quantity, and packed close together, in which case it can be treated as a continuous solid. Treating the lung parenchyma as a homogeneous elastic solid has been previously done by Kamm [14], Lai-Fook [15] and Ma [27]. It can be seen in figure 2.1 that this assumption is valid, since bronchiole is encased within the lung parenchyma tissue.
However, on a small scale, the lung parenchyma seems to be more appropriately described as a discrete spring network, as seen in figure 2.5. In the second solution of the model, the lung parenchyma will be modeled as a discrete network of springs that surrounds the airway. This network of springs, which represents the lung parenchyma, resists contractile deformation and prevents full closure of the airway. In this case, the lung parenchyma will be studied at a certain fixed volume, and previously obtained experimental results for the mechanical properties of the lung parenchyma will be used to determine the effects on airway contraction.
It is the aim of this project to obtain values for the percent change in airway caliber, as well as to understand the deformation throughout the airway. The percent change in the airway caliber is hypothesized to depend on the geometrical dimensions of the airway, as well as on the inner pressure of the airway, and finally the elastic properties of the layers. To obtain meaningful results from the models developed, general values of the physical geometry must be known.  In physiology, airway generation means a division point in the airway in which one airway branches into two or smaller airways. Generation 0 is the trachea, or windpipe. Generation 16 is the terminal bronchioles. It can be seen in table 1 that there is a jump in diameter between generation 4 and 5, and then a steep drop in diameter once generation 6 is reached. This is to be expected, since the human airway generation will have different values of geometry, depending on the lung in question.
According to Codd [4], the values of Young's modulus for the smooth muscle layer averaged 3.3 kPa for circumferential strains. Codd also estimated the Young's modulus for the stiff layer to be 20 kPa in the circumferential direction. The shear modulus in terms of Young's modulus and Poisson's ratio is The Poisson's ratio in this case is 0.5 due to the incompressibility of the materials. The shear moduli for the outer and the stiff layer can be easily calculated to be 1.1 kPa and 6.67 kPa, respectively. The shear modulus for the lung parenchyma is given by Where is the transmural pressure. The transmural pressure will be taken at the functional residual capacity, and is usually given to be 0.49 kPa [21]. This gives the shear modulus for the lung parenchyma to be 0.343 kPa. Figure 2.7 shows the airway layers with their respective material properties labeled.

SECTION 3 -DEFORMATION
Since the geometry is cylindrical, it is convenient to develop the analytical model in cylindrical polar coordinates ( , , ). Consider the two layer cylinder to be stress free, as the initial or reference configuration. The geometry can be described as Where is the inner radius of the stiff layer, B is the outer radius of the stiff layer, and is the outer radius of the smooth muscle layer, as seen in figure 3.1. L is the length of the airway. Note that in figure 3.1, the schematic shown is that for the lung parenchyma modeled as an infinitely large domain. Since the mass in this system is conserved, it implies that every particle in the reference configuration has a one to one mapping to the current configuration. This deformation mapping function maps the reference position of every particle to its current position . The deformation mapping function is extended and can be expressed as the deformation tensor, denoted by Where the operator ∇ 0 is simply the gradient Taking the gradient of the position vector in the current configuration gives Where the subscripts , and are the partial derivatives , and , respectively.
For an axis-symmetric long tube, the deformation gradient reduces to The deformation gradient maps a material vector in the reference configuration, to a spatial vector in the current configuration, and vice versa, of the neighborhood of a material particle. The deformation gradient includes change in lengths, or stretching, as well as changes in angles, or shearing. However, the deformation gradient also includes a rotation of the neighboring particles. Since rotation does not play a role in shape change, it is useful to decompose into two tensors. One that accounts for shape change, and the other for rotation. Imagine two scenarios of the deformation gradient. One in which the deformation gradient only describes pure rigid body displacements, without any deformations. Then the deformation gradient can be represented by a proper orthogonal tensor . The other deformation gradient describes pure stretching and in this case is a symmetric positive definite tensor , or . Finally, the deformation gradient can be expressed by The equation above describes the polar decomposition of the deformation gradient. In equation (3.7), is the right stretch tensor, and is the left stretch tensor. Taking equation (3.7) further, it can be said that The deformation gradient gives a one-to-one mapping, for the position of a neighborhood of particles in the reference configuration, to the current configuration, . This one-to-one mapping imposes a restriction on the deformation, that is, particles are neither created nor destroyed. It means that a single particle cannot be mapped to two positions, and that two particles cannot be mapped to one position. Mathematically this can be said that ≠ . Since > 0, then the term in equation (3.3) cannot be equal to zero, unless the term is zero. Therefore, • is a positive definite quadratic in the form of: Where is a symmetric positive definite tensor. This means that its square roots exist. is the right Cauchy-Green deformation tensor. From equation (3.9) it can be seen that the right Cauchy-Green deformation tensor can be mathematically described Using equation (3.6) the right Cauchy-Green deformation tensor can be expanded For the special case of axis symmetry and plain strain assumption, it is assumed that the displacement is given by = ( ). this means that the deformation and displacement is independent of , and =0, therefore , and , are zero.
The deformation gradient , and the right Cauchy-Green deformation tensor reduce to Since the right Cauchy-Green deformation tensor is symmetric, its eigenvectors are automatically mutually perpendicular as long as no two eigenvalues are the same. If two or all three eigenvalues are the same, the eigenvectors are not uniquely defined. In this case, any convenient mutually perpendicular set of eigenvectors can be used. It follows immediately from equation (3.12) that the principal axes of are , , .
can then be expressed in terms of eigenvectors and eigenvalues The square root of gives the stretches in the principal directions, that is It is assumed that the materials are incompressible so that the volume is preserved.
Here is the Jacobian. The definition of the Jacobian leads to a local condition for invertibility. As previously mentioned, the deformation gradient is a one-to-one mapping in the neighborhood of the particle .

SECTION 4 -HYPERELASTICITY
In the paper by Codd [4], mentioned earlier, the strips of membranous trachea were strained to roughly 20%. Biological materials can obtain strains much larger than that, it can be from 50-100% strain. The first Piola Kirchhoff stress tensor and the Cauchy stress tensor can be written as The invariants of the right Cauchy-Green deformation Tensor are a function of the right stretch tensor, that is 1 ( 2 ), 2 ( 2 ), where is the right stretch tensor.

-INFINITE DOMAIN
It is possible for the airway to experience an internal on the periphery of the stiff layer.
Let 1 be the inside pressure, acting on the stiff layer. The boundary condition for the elastic tube under external pressure is given by The lung parenchyma layer is also modeled as an infinite solid continuum. In such a case, equation ( Where is the stress exerted by the lung parenchyma on the airway. In this case, there will be three integration constants to solve, one for each individual layer.

-SPRING DOMAIN
As previously mentioned, the parenchymal attachments act as springs, with spring constant ̅ , which prevent deformation of the airway. The second boundary condition aims to incorporate the effect that the lung parenchyma has on the caliber change. The stress on the boundary of the smooth muscle layer exhibits the stress that the lung parenchyma exerts, however in this case; it is modeled as a series of linear springs that undergo displacement. This displacement changes the magnitude of the stress exerted by the lung parenchyma, which is exactly what is seen in an actual lung, since the stress exerted by the lung parenchyma is a function of the lung volume.

-SOLUTION FOR THE STIFF LAYER
As shown in equation (   This solution gives the Cauchy stress as a function of the radius in the current configuration. Of course, it should be noted that this solution is not completed because the deformed inner radius is still unknown.

-SOLUTION FOR THE SMOOTH MUSCLE LAYER
Upon being stimulated, the smooth muscle layer of the airway contracts and The active stress in this case is developed in the direction, therefore the Cauchy stress component in the direction will be broken down into an active and passive term. The total Cauchy stress component (2) is given by Where the subscripts and are the active and passive terms, respectively. As previously mentioned, all things with the superscript (2) represent all things relating to layer two. The passive term is already given by equation (5.4). Substituting this value, the total Cauchy stress in the direction, is given by  Note that the solution for layer two was not included in equation (6.12). This is solely for simplicity. Finally, substituting the value for (2) and again using the relationship between reference and current radii gives 2 )] + (2) ( ) − 1 + ̅ = 0 (6.12) The above equation fully relates the active stress to the displacements in the airway.
The desired results are for the deformed inner radius to be a function of the active contractile stress (2) that is = ( (2) ). The radii in the reference configuration , , , are assumed to be known. The radii in the current configuration and are directly related to the inner radius by equations (3.18). Therefore, equation (6.13) gives the contractile stress as a function of the inner radius in the current configuration. Equation is a nonlinear, logarithmic algebraic equation, making it difficult to obtain a solution for in terms of (2) . To solve the above equation in order to obtain the caliber change, a MATLAB program has been completed, employing the bisection method.

-INFINITE DOMAIN SOLUTION
In the solution for the smooth muscle layer, the term ̅ represents the effect of the This solution takes into consideration the stress developed by the lung parenchyma in response to the displacement in the airway. As previously mentioned, the mechanical properties of the lung parenchyma are dependent on the transmural pressure [15].

SECTION 7 -THE BISECTION METHOD
Certain solutions can be obtained with ease without the need of numerical analysis. For example, if is a polynomial of degree 4 or less, formulas for its roots exist. However, if the degree of the formula is 5 or higher, with arbitrary coefficients, then it does not have a general algebraic solution, as stated by the Abel-Ruffini theorem [16]. Also, if is a general nonlinear function, then no closed-form algebraic solution exists, as is the case with equation (6.13). This entails that a method that computes approximate roots must be utilized in order to obtain a solution for this equation. Considering a function ( ) = 0. An iterative method can be used to compute its solution numerically, by construction a sequence of 1 , 2 , … , , +1 , … that converge to a root of ( ) = 0 . There are several mathematical methods used in order to converge this sequence to the root of ( ) = 0, such as the Secant method, Newton's method, or bisection method. The method used here is that the bisection method.
A bisection is the division of a given curve, or an interval, into two equal parts. , or bisect the interval, to subintervals that contain the value 0 . This process can be performed several times, locating the half of the interval that contains 0 . Note that the procedure will work even if there is more than one root in the interval [ , ]. However, for simplicity it is assumed that there is one unique root. The bisection method is slow and time consuming. For the purpose of this project however, a faster method is not required since the solution necessary to produce satisfactory results does not require a large number of iterations.
In more detail, a bisection procedure is used in iterations, which converses on a solution that is known to lie inside of the given interval, [ , ]. The solution to the function is given by To start, set the first interval as 1 = and 1 = and let 1 be the midpoint of [ , ].
For the first iteration, the solution is Note that if ( 1 ) = 0 then no further bisection is necessary, and the solution has been located. If ( 1 ) ≠ 0 then ( 1 ) has either a positive or negative value. If ( 1 ) and ( 1 ) have the same sign, then the solution is found in the interval [ ]. Then set 2 = 1 and 2 = 1 . In either case, the interval is half as long as the initial interval. The process selects a subinterval that guarantees to have a root. The next interval will then be [ 2 , 2 ]. The advantage of the bisection method is that it is very straight forward and simple, as well as it always provides the roots of a given continuous function with the property ( ) ( ) < 0.
If 1 is the midpoint of the initial interval, and is the midpoint of the ℎ interval, then the difference between and a solution is bounded by Where is a given error, or tolerance. Taking the natural logarithm of both sides gives Equation (7.4) can be used to determine in advance the number of iterations that the bisection method needs in order to converge to a root within this tolerance [17]. A pseudocode for the bisection method is given in Appendix II, as well as the MATLAB code that solves for the analytic solutions.

-LINEAR BUCKLING
In a study by Klingele & Staub [18], the terminal bronchioles in cat lungs were studied, in which they freeze dried the lungs after inflating the lungs to 100%, 48%, 23%, and 8% of their total lung volume. The diameters were then measured. They observed that at the lowest lung volume, the airway walls had buckled into kinks and folds. Since the buckling of airway wall is a well-known and documented phenomenon, it is necessary to obtain reliable results when the buckling is included.
Buckling is a mode of failure that is due to elastic instability of an equilibrium configuration, without fracture of separation of the continuum. Buckling, often referred to as loss of stability, refers to a large increase in displacement when there is a small increment in load. In most structures, the displacements increase gradually with increased applied load. If the applied load is too large, especially for compressive structures, a small increase in applied load can lead to a sudden large increase in the displacements. Buckling refers to this transition to large, often catastrophic displacements. Buckling mainly occurs due to thermal or mechanical loads. Buckling is a type of instability. Instability is the state in which small perturbations can cause large changes in the response of the structure. The instability in buckling is due to geometric effect and not a change in the material property, as is in material yielding and failure. In this course we will limit ourselves to instability of structures due to buckling. The differential equation for the loss of stability in columns is Where is the modulus of elasticity, is the moment of inertia, is the non-trivial equilibrium mode of the buckled column, and 0 is the axial compressive load [19].
This differential equation represents an eigenvalue problem.
An eigenvalue problem is defined as one in which the values of the parameter are obtained from the equation The obvious solution for equation ( The ABAQUS theory guide makes use of several continuum mechanics identities and constitutive equations in order to rewrite the equation into a form that can fully describe elasticity, hypoelasticity, and hyperelasticity.
Using the finite element approach, the governing equation for buckling presented above takes the form of Where 0 is the stiffness matrix in the reference configuration, which includes the effects of the preloads, ∆ is the differential initial stress and load stiffness matrix, are the eigenvalues, and are the eigenvectors, the non-trivial displacement solutions, often referred as the buckling mode shapes. The magnitude of the eigenvectors is normalized in a linear buckling analysis. They do not represent the magnitude of deformation at the critical load. The buckling mode shapes do not have a magnitude since they are used to predict the likely failure mode of the structure. For this fully formulated finite element solution, the applied loads can consist of pressures, concentrated forces, nonzero prescribed displacements, and thermal loading. The loading in which the system will buckle, called the critical buckling load, is then Where is simply the result from the applied forces and , as well as the prescribed displacements , and is the perturbation loads, due to ∆ ∆ and ∆ .

-POST BUCKLING
The problem described is that of a nonlinear, hyperelastic multilayered cylindrical airway. The buckling analysis performed by ABAQUS is that for an elastic buckling by eigenvalue extraction. The estimation is useful for stiff structures.
However, the problem formulated is nonlinear, and cannot be fully described by the linear buckling analysis. In the case in which the structure responds in a nonlinear way prior to buckling, as is the case, a nonlinear analysis using the Riks method is required in order to obtain reliable results. The Riks method is used to predict the nonlinear collapse of a structure. The Riks method solves for loads and displacements, simultaneously. In order to do so, ABAQUS uses an "arc length" , in a loaddisplacement space, in order to measure the progress of the solution. The Riks method in ABAQUS uses the solution for the buckling modes from the linear buckling analysis. Since this is a continuation of a previous history, the loads from the previous analysis are not redefined and are treated as "dead" loads 0 , with constant magnitudes. The load defined during the Riks step is referred to as a "reference" load, . Since the load varies through the load-displacement space, the current load must always be proportional to the dead and reference load. The total load is given by The value should not be confused with the eigenvalue given previously. Here is the load proportionality factor. This load proportionality factor is found as part of the solution given by the Riks method. The load proportionality factor is obtained at each increment in the load-displacement space. `The initial load proportionality factor, ∆ is computed as Where is a user specified total arc length [20].

-INFINITE DOMAIN
For generation 0, the radii to the inside and outside of the stiff layer are 8410.6 μm and 8578.3 μm, respectively.  As can be seen in figure 9.1, the trend is similar to the general graph of a negative natural logarithm function. As expected, figure 9.2 follows a very similar trend, showing similarities to the exponential graph of .  25.51 percent, justifying the large deformation solution. It can also be seen that the difference in percent caliber change between generations 0 and 16 is 1.65 percent.

-SPRING DOMAIN
In order to obtain results for the spring domain, a value for the spring constants , ) The Buckingham Pi theorem states that the total number of these relevant dimensional parameters (n) can be grouped into n-m independent dimensionless groups. The number m is usually equal to the minimum of independent dimensions required to specify the dimensions of all relevant parameters. Since (1) and (2) have the same dimensions as (2) , these parameters will not be used in order to obtain the dimensionless parameter. The dimensions here are ̅ = 3 • 2 (2) =

SECTION 10 -FINITE ELEMENT SIMULATION
As was shown in the previous sections, the generation of the airway has a buckling critical load. Once this load is reached, the structure will buckle, and the axisymmetric deformation assumption will no longer be valid. These values will be considered when performing the numerical analysis.
In order to simplify the simulation for finite element analysis, only the infinite domain will be modeled due to the fact that it describes the problem more accurately than the spring domain. Generation 0 (trachea) is modeled using ABAQUS and presented here. To successfully complete the ABAQUS simulation, first a linear buckling analysis must be completed. The linear buckling analysis is necessary since it predicts the buckling mode of the structure, and gives the critical buckling load that will cause the structure to buckle. The mathematical model developed assumes plane strain, axisymmetric deformations. These assumptions carry to the finite element simulation as well.  The loading will be simulated as a thermal contraction. Figure   As can be seen, the soft muscle layer is the only layer which has a thermal load applied. The high temperature is given to be 0, built into the initial step of the simulation. The low temperature is given to be -1, built into the Buckle step of the simulation. It is only necessary to provide a temperature difference of magnitude 1, since the critical buckling load is given by the eigenvalue multiplied by the prescribed load, as seen in equation (8.5). The coefficient of thermal expansion for the material of layer two is 1 in order to simplify the results.
A boundary condition is applied on the periphery of the cylinder. A cylindrical coordinate system is created in order to fix the and directions. Since the third layer is assumed to be infinite, then it will have zero stress, strains and displacements on the far boundary, and the infinite domain is simulated here approximately by a large domain and the outmost surface are held fixed for the simulation.   Different mesh sizes were attempted prior to concluding that the size of 0.0001 provided the desired results. As can be seen in figure 9.5, the mesh is evenly constant throughout the geometry.
Before submitting the job for completion, it is important to create an .fil file. The .fil file contains the geometry file of the buckling analysis. This file contains the imperfections in the geometry from the buckling. In order to create a .fil file, right clicking the "Part" tab on the left of the screen, the option for "Edit Keywords" is selected. At the end of the keyword file, the command *NODE FILE U must be written, as seen in figure 10.5. Once this is completed, the job is submitted with the file name thermal_buckle.  This command obtains the imperfections from the .fil file created. Where 1 in the command given is the buckling mode 1, which will be used for the imperfection, and 1e-4 is the degree of the imperfection to be included.        enough, and the trend could be a logarithmic trend. This implies that the critical contractile buckling stress will increase with increase in the overall geometrical dimensions of the cylinder. It is also worth noting that the structure will buckle according to the material properties given to each layer.  Similarly, the inner pressure term can be removed and the solution will still hold true.
Moreover, it is quite straightforward to only complete the solution for a single layer if it is desired. This formulation covers a several solutions for cylindrical problems. The formulation of the model can be applied for the lumen of the digestive track. Figure 11.1: Schematic illustrating the form of ridges and folds in the lumen of the digestive track. Figure reproduced from [26] As seen in figure 11.1, the digestive track also implements smooth muscle contraction. The solution presented in this paper can be used in studies for airway tissue engineering. An in vitro model based on the human airway musculature-on-chip technology has been developed [23]. The human airway musculature on chip technology consists of human bronchial smooth muscular thin films, comprised of a bottom layer of elastic polymer polydimethylsiloxane (PDMS), and a top layer of engineered bronchial smooth muscle. It is capable of simulating healthy and asthmatic bronchoconstriction and bronchodilation. When stimulated, the smooth muscle layer contracts and bends into a curvature, which is analogue for the constriction of the human airway.  As can be seen in figure 11.2, the human airway on a chip is a flat cantilever, and does not include a real world representation of the caliber change in the airway.
In order to extrapolate the results of the musculature-on-chip flat cantilever to a real world representation of the caliber change, the theoretical model presented in this document can be used.
The finite element simulations give a better understanding of the behavior which the airway exhibits when it is under contractile stress. The simulations provided values of contractile stress for which the theoretical model will fail. The simulations presented here seem to be a better means of simulating airway contraction, since the contractile loads applied are more closely related to those seen in real world physiology. Also, the modeling of the lung parenchyma as an infinite domain is able to capture their effect on the airway deformation, and provide good results of the critical contractile stresses for the use in the theoretical model.

SECTION 12 -FURTHER WORK
To fully complete this work, if the assumption that the lung parenchyma is modeled as an infinite domain will be carried on to the finite element analysis, then the elements used to model the lung parenchyma will need to be modeled as infinite elements, and not finite elements as done in this work. It is also necessary to write a code for ABAQUS in order to incorporate a constitutive model that can more accurately describe the contractile stress, instead of modeling the contractile stress as thermal contraction. To obtain a more accurate mathematical model, it is possible to drop the assumption of isotropy, and model the materials as anisotropic, since this seems to be a more accurate way to model biological materials. Also, it should be mentioned that during abnormal conditions, the airway also swells and grows. This However, the most important portion that is necessary in order to complete this work is to have experimental results to validate the mathematical model. Due to the nature of the problem, obtaining reliable results for the caliber change after a certain prescribed stress has been introduced is extremely difficult. In a physiological organism, knowing the exact contractile stress generated by the smooth muscle is difficult to obtain, as well as measuring the caliber once the stress has been induced.
However, this is necessary in order to validate the model. Step 1 i = 1; FA = f (a).
Step 2 While i < N0 limit iterations to prevent infinite loop. Do Steps 3-6.