Abstract
Pulmonary diseases, driven by pollution, industrial farming, vaping, and the infamous COVID-19 pandemic, lead morbidity and mortality rates worldwide. Computational biomechanical models can enhance predictive capabilities to understand fundamental lung physiology; however, such investigations are hindered by the lung’s complex and hierarchical structure, and the lack of mechanical experiments linking the load-bearing organ-level response to local behaviors. In this study we address these impedances by introducing a novel reduced-order surface model of the lung, combining the response of the intricate bronchial network, parenchymal tissue, and visceral pleura. The inverse finite element analysis (IFEA) framework is developed using 3-D digital image correlation (DIC) from experimentally measured non-contact strains and displacements from an ex-vivo porcine lung specimen for the first time. A custom-designed inflation device is employed to uniquely correlate the multiscale classical pressure-volume bulk breathing measures to local-level deformation topologies and principal expansion directions. Optimal material parameters are found by minimizing the error between experimental and simulation-based lung surface displacement values, using both classes of gradient-based and gradient-free optimization algorithms and by developing an adjoint formulation for efficiency. The heterogeneous and anisotropic characteristics of pulmonary breathing are represented using various hyperelastic continuum formulations to divulge compound material parameters and evaluate the best performing model. While accounting for tissue anisotropy with fibers assumed along medial-lateral direction did not benefit model calibration, allowing for regional material heterogeneity enabled accurate reconstruction of lung deformations when compared to the homogeneous model. The proof-of-concept framework established here can be readily applied to investigate the impact of assorted organ-level ventilation strategies on local pulmonary force and strain distributions, and to further explore how diseased states may alter the load-bearing material behavior of the lung. In the age of a respiratory pandemic, advancing our understanding of lung biomechanics is more pressing than ever before.
Introduction
Respiratory diseases and disorders, such as asthma, emphysema, bronchitis, pulmonary fibrosis, and lung cancer, collectively lead as the global cause of morbidity and mortality (Centers for Disease Control and Prevention, 2015; Eskandari et al., 2018). These pulmonary illnesses impose strenuous social and economic burdens, as seen with the recent lung-damaging COVID-19 outbreak (Atkeson, 2020). The acute and progressive pathological inflammation and bronchoconstriction of the lung obstruct and restrict airflow and oxygenation, inducing altered mechanical properties (Suki and Bates, 2011; Eskandari et al., 2016). This mechanical remodeling is multiscale, spanning the destruction of alveolar sac elasticity in emphysema (Suki et al., 2003), the over stiffening of the parenchymal tissue in pulmonary fibrosis (Faffe and Zin, 2009), and the constriction and collapse of airways in asthma (Bai and Knight, 2005; Eskandari et al., 2015; Maghsoudi-Ganjeh et al., 2021). Thus, the hierarchical and complex structure of the lung highlights the importance of mechanics in respiratory health (Tawhai and Bates, 2011; Eskandari et al., 2013).
Despite the growing body of literature on pulmonary mechanics, the multiscale and multiphysics link between the global pressure-volume behavior of the lung and the local-level tissue deformation remains largely unexplored. There has been notable progress to characterize the lung at the organ scale through classical pressure-volume curves and at the tissue level using indentation and uniaxial tensile tests (Lai-Fook et al., 1976; Zeng et al., 1987; Fung, 1988; Eskandari et al., 2018); however these investigations remain siloed at disconnected scales. Amalgamating these multiphysics and multiscale behaviors is central to understanding lung disease mechanisms, predicting disease progression, and mitigating ventilator-induced-lung-injuries (VILI) to eliminate tissue over stretching (volutrauma) and stressing (barotrauma) (Dreyfuss and Saumon, 1998; Vlahakis et al., 1999; Arora et al., 2017; Arora et al., 2021). Unless an atlas for pulmonary kinetics and kinematics can be established, current ventilation protocols will continue to be subject to trial and error approaches and hindered from advancements despite exigent demands instilled by a worldwide pandemic (The Acute Respiratory Distress Syndrome Network, 2000; Amato et al., 2015).
Advancements in biologically-oriented digital image correlation (DIC) techniques have facilitated quantifying the mechanical connections between organ-level breathing and local tissue behavior for fast, large, and non-linear deformations. DIC is a common full-field, non-contact deformation characterization technique originally applied on inert structures (Chu et al., 1985), and has now been enhanced to study the behavior of intricate biological tissues, such as the cornea (Boyce et al., 2008), arteries (Sutton et al., 2008), knees (Mallett and Arruda, 2017), and most recently, the lung (Mariano et al., 2020). In this method, sequential images of a specimen’s speckled surface undergoing loading are used to obtain the topological displacement field (Chu et al., 1985). While DIC describes the kinematics, inverse finite element analysis (IFEA) can be employed to divulge the kinetics. IFEA yields specimen mechanical properties by minimizing the error between the displacements predicted by the Finite Element (FE) model and those measured via experiment (Birzle et al., 2019).
Here we construct the first in-silico IFEA structural representation of the whole lung as informed and validated from DIC resulting from applied evolutionary pressure-volume loading controlled by a custom-designed breathing apparatus (Mariano et al., 2020; Sattari et al., 2020). Based on the obtained surface geometry and deformation map of the inflating lung, a corresponding reduced-order 3-D FE model is constructed using membrane elements undergoing the same experimental lung pressures. Various constitutive models are explored, including homogeneous isotropic hyperelastic, homogeneous anisotropic hyperelastic, and heterogeneous isotropic linear-elastic materials (Mooney, 1940; Holzapfel et al., 2000). The parameters of the multiple FE models are calibrated through a fully automated IFEA framework. Both classes of derivative-based and gradient-free optimization algorithms are implemented to predict the material response by minimizing the error between model-predicted and DIC-recorded displacements of the external surface of the lung. The set of calibrated material parameters, along with local and heterogeneous deformation results of the in-silico lung, are presented, model performances are compared, and future applications are discussed.
Materials and Methods
Digital Image Correlation and Pressure-Volume Experiments
Previously established extensive experimental DIC protocols and pressure-volume tests were utilized for the ex-vivo specimen tests conducted here and will be briefly summarized (Mariano et al., 2020; Sattari et al., 2020). Fresh porcine lungs from an abattoir were obtained (50 kg female domestic York farm minipigs, Institutional Animal Care and Use Committee approval not required) and a plastic tube was inserted through the trachea to fully inflate the lung using an airline pressure system. A generic exfoliator pad dipped in quick-drying white enamel paint (rust-oleum) was used to create speckles (Mariano et al., 2020). The specimen was loaded into our custom pressure-volume apparatus for controlled inflation tests; this device consisted of two pistons (a source and a response), a transparent tank, and a computerized controller system (Sattari et al., 2020). 900 ml of air was applied to the lung, and the real-time continuous pressure and actual volumetric deformation of the lung (less than the applied volume due to the compression of air) was measured. As in previous studies, a preload of 5 cmH2O was used as the reference state. A rate of 15 breaths-per-minute was used and the specimen was preconditioned three times to generate reproducible cycles and the fourth inflation response was analyzed.
The 3-D stereoscopic DIC system (ARAMIS 12M, Trilion Quality Systems) consisted of two optical cameras hovering over the transparent tank, which recorded the dynamic deformation of the lung at 10 Hz. For a measuring volume of 375 × 295 × 295 mm, the calculated displacement measurement accuracy was 0.10875 mm (Jones and Iadicola, 2018). The images were analyzed following standard DIC techniques to calculate the displacement and strain of the exterior surface of the lung relative to its uninflated state (Chu et al., 1985; Jones and Iadicola, 2018) and are extensively detailed in Mariano et al., 2020. Figure 1A showed the DIC strain map corresponding to the peak pressure-volume inflation stage. Figure 1B showed the corresponding pressure-volume curve obtained from the inflation. Based on this curve and the inflation rate, the pressure-time amplitude curve (Figure 1C) was extracted and applied to the FE model as the loading step.
FIGURE 1
Inverse Finite Element Analysis Overview
The FE models were calibrated using the measured DIC displacement of the ex-vivo lung. The algorithm (shown in Figure 2) minimized the error by finding the optimal values for unknown free material parameters. The inverse FE model leveraged several known parameters: the surface geometry of the lung from the uninflated stage obtained after three preconditioning inflation-deflation cycles, the experimentally measured pressure-time graph, and the DIC displacement field of the lung surface. However, the type of constitutive model and corresponding material parameters were unknown; these model attributes were left to be determined by the optimization algorithms. Before applying the model to the actual lung, we verified that when applied on a simpler geometry with known deformation field and material properties the model is able to recover the given material parameters successfully.
FIGURE 2
Various constitutive relations were examined. The material parameters were initialized from several distinct starting points and the solution was generated multiple times to ensure mathematical robustness. The error was calculated based on the normalized squared sum of residuals. Material properties were perturbed to calculate the sensitivity of the error, which informed the alternate directions adopted by the optimization algorithm. This incremental procedure was repeated to progressively minimize the error until a pre-defined convergence criterion was met (change in the error less than 10−6 of the initial value). Given the nonlinear lung pressure profile (Figure 1C), the displacement error was evaluated at several increments evenly spaced out throughout the inflation cycle and not just at the full inflation point.
Finite Element Model Organization
The built-in stereo camera DIC system was used to capture the exterior surface of the lung. The obtained geometry was tessellated with 3-D triangular elements averaging 1 mm in size. The DIC system could only detect and analyze the visible portion of the lung lobes (Figure 1). The raw data from the original DIC geometry was not suitable for the FE simulation as it contained some elements with poor isoperimetric quality and sharp surface discontinuities. To improve the mesh quality, MeshLab (ISTI-CNR, Italy) was used to smoothen the surface while preserving the original surface features (Cignoni et al., 2008). Two fine and coarse meshed models, with ∼5,000 and 457 elements respectively, were exported as STL files which were then converted to Abaqus input files using the built-in script plugin (Dassault Systems, Providence, RI, United States). The fine and coarse meshes were used to study the cases corresponding to homogeneous and heterogeneous material models, respectively. This approach reduced the number of unknown material parameters in the heterogeneous case and was necessary to substantially decrease the IFEA computational cost. After confirming the mesh resolution was sufficient, both meshed models were discretized using 3-D membrane elements (M3D3). The recorded DIC experimental displacement values of the nodes sitting on the perimeter of the surface geometry were applied as displacement boundary conditions to the FE models; as such, while the geometry of the model represented the visible portions of the lung lobes, the role of the adjacent tissue was represented through the application of these periphery nodal displacements (as opposed to traditional boundary conditions). The thickness of membrane elements were set to 1.0 mm.
One-to-one experimental to FE model nodal correspondence was created for error calculations by probing the displacement from ∼7,000 evenly distributed points across the surface. During the multiple inflation stages, where the displacement error between the FE model and DIC nodes were to be calculated, there is not a one-to-one correspondence between FE model nodes and DIC probe points. Therefore, interpolation must be used to find FE model nodes corresponding to the DIC probe points so that the error could be computed. An interpolation technique was used, utilizing the DIC displacements and coordinates of the probe points to train the k-nearest-neighbor algorithm (Altman, 1992). The performance of the interpolation model was evaluated using 10-fold cross validation, confirming the accuracy was above 0.95. The number of nearest neighbors k was set to be five. The interpolation technique was implemented using two Python scripts: one to access nodal coordinates from Abaqus, and one to perform interpolation.
The FE model was solved using Abaqus dynamic implicit solution scheme by subjecting the model to the experimental lung pressure values at five deformation stages (Figure 1C). Nonlinear geometry formulation was utilized, and the displacement and strains were analyzed. In order to significantly accelerate the IFEA process, all the scripts were parallel-coded such that multiple FE simulations were running simultaneously in batch mode.
Constitutive Models
The best performing constitutive model was not known a priori. We investigated three different material model cases to consider homogeneity versus regional heterogeneity, preferential orientation using an anisotropic versus isotropic response, and the linear versus nonlinear cases to determine optimal constitutive parameters as detailed below. It was important to note the reduced-order nature of the model meant the parameters of these constitutive models were not simply the material properties of the lung; rather they are a pseudo-material model, represented as a projected, averaged surface response of tissue or compound material parameters of the parenchyma, airways, and pleura layer consolidated together.
Homogeneous Isotropic Hyperelastic Case
Here the compressible Mooney-Rivlin hyperelastic model (Mooney, 1940) with the strain energy density defined as was used. , were the first and second invariants of the deviatoric deformation tensor, was the Jacobian of deformation gradient , and and were the three unknown material parameters. The stress-stretch curve of this model can span strain-hardening and strain-softening behaviors, depending on the relative values of the material parameters, allowing a versatile IFEA framework.
Homogenous Anisotropic Hyperelastic Case
The Holzapfel-Gasser-Ogden (HGO) formulation (Holzapfel et al., 2000) with the strain energy density defined as was used here. and were defined as = and . The five unknown material parameters were and . The three strain-representing kinematic variables were , , and . In this formulation was the first invariant of the deviatoric deformation tensor. In addition, was the pseudo-invariant of and where followed the multiplicative decomposition of the deformation gradient and the deformation tensor . The vector was a unit vector field defining the fiber direction in the undeformed configuration. was the Jacobian of deformation gradient . represented the squared of the stretch ratio of the material fiber in the direction of the fiber family defined by . The degree of preferential alignment of the fiber family governing anisotropy was controlled by the dispersion parameter , where 0 indicated the family of fibers were fully aligned and 1/3 indicated a completely random distribution of the fiber family (reducing to isotropic form). For each element in the mesh, the local z-direction was the outward normal to the element surface, the local y-direction was specified in the anterior-posterior direction, and the local x-direction was subsequently determined based on the right-hand sign convention for Cartesian coordinates. Given that strains were smaller in the medial-lateral direction, was aligned with the defined local x-direction for each element. It should be noted that even though we define x-axis for the main fiber directions, since the parameter is left free to be determined by the optimization engine, the true fiber directions are not strictly fixed in the model.
Heterogeneous Isotropic Linear-Elastic Case
In this case we considered a linear elastic isotropic material model (Eskandari and Kuhl, 2015) with Young’s modulus and Poisson’s ratio . The regional heterogeneity of this model meant each element of the mesh had its own two parameters that were found by the optimization algorithm. The shear modulus of each element was then calculated using and regionally mapped onto the lung surface.
Optimization Algorithms for Extracting Material Parameters
Gradient-based optimization algorithms are prone to returning local optima in the neighborhood of the initial search point, while the derivative-free optimization algorithms, such as meta-heuristic algorithms, are more likely to return the global optimal solution instead (Nocedal and Stephen, 2006). Therefore, two broad classes of gradient-based and gradient-free optimization algorithms were implemented to improve global optimum acquisition. In the gradient-based approach, the trust-region-reflective (TRR) algorithm (Steihaug, 1983), available in Matlab lsqnonlin function (The MathWorks Inc., MA, United States), was used for the two homo/iso/hyper and homo/aniso/hyper cases; and the sequential-quadratic-programming (SQP) algorithm (Boggs and Tolle, 1995), available in Matlab fmincon function, was used for the hetero/iso/linear-elastic case. As for the gradient-free approach, we implemented our own version of the particle swarm optimization (PSO) algorithm (Kennedy and Eberhart, 1995) in a Matlab script. This provided more flexibility to impose custom constraints on our problem, such as bounds on the position and velocity of particles, which may not be done so freely in the built-in Matlab PSO algorithm. The IFEA processes in both approaches were fully automated by conjoining Matlab, Abaqus, and Python. To help avoid local optima, the TRR and SQP algorithms were run from several randomly selected initial estimate points within the search space as given in Table 1 and Table 2. As a verification step, we applied the IFEA framework to a test case model with a simpler geometry and known deformation field and material parameters. We confirmed that the optimization pipeline was indeed successful in recovering the pre-known material parameters.
TABLE 1
| Initial material parameters | Optimal material parameters | |||||
|---|---|---|---|---|---|---|
| Run # | (kPa) [1–200] | (kPa) [1–200] | [10−4–10−2] | (kPa) | (kPa) | |
| 1 | 25 | 10 | 5 | 136.2 | 1.0 | 13.2 |
| 2 | 10 | 20 | 15 | 136.6 | 1.0 | 13.5 |
| 3 | 100 | 42 | 61 | 136.6 | 1.0 | 13.5 |
| 4 | 172 | 18 | 32 | 136.9 | 1.0 | 13.6 |
| 5 | 10 | 90 | 80 | 136.1 | 1.0 | 13.2 |
| 6 | 50 | 50 | 50 | 136.7 | 1.0 | 13.6 |
| 7 | 200 | 150 | 90 | 136.5 | 1.0 | 13.4 |
Sets of initial estimates and converged optimal material parameters for the homo/iso/hyper case. In order to avoid local minimum, the optimization routine was repeated from seven different starting points.
TABLE 2
| Initial material parameters | Optimal material parameters | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Run # | (kPa) [1–200] | (×10−4)[10−4–10−2] | (kPa) [1–200] | [0–1] | [0–0.33] | (×10−4) | ||||
| 1 | 18 | 42 | 8 | 0.30 | 0.11 | 116.5 | 42 | 1.0 | 0.05 | 0.33 |
| 2 | 51 | 28 | 14 | 0.70 | 0.05 | 116.6 | 28 | 1.0 | 0.24 | 0.33 |
| 3 | 80 | 12 | 100 | 0.45 | 0.2 | 116.2 | 12 | 1.0 | 0.09 | 0.33 |
| 4 | 150 | 10 | 76 | 0.62 | 0.3 | 116.0 | 10 | 1.0 | 0.13 | 0.33 |
| 5 | 7 | 100 | 35 | 0.90 | 0.17 | 116.3 | 100 | 1.0 | 0.15 | 0.33 |
| 6 | 180 | 63 | 22 | 0.10 | 0.17 | 116.5 | 63 | 1.0 | 0.07 | 0.33 |
| 7 | 100 | 50 | 150 | 0.18 | 0.25 | 116.5 | 50 | 1.0 | 0.13 | 0.33 |
Sets of initial estimates and converged optimal material parameters for the homo/aniso/hyper case. In order to avoid local minimum, the optimization routine was repeated from seven different starting points.
The Adjoint Method to Calculate the Objective Functional Gradient
In this study, the adjoint method was used to calculate the gradient for the hetero/iso/linear-elastic case given that the total number of unknown material parameters for this case was much greater than that of the homo/iso/hyper (three parameters) and homo/aniso/hyper (five parameters) cases; In the former case there were 457 elements and two material unknowns at each element, which rendered the total number of unknown parameters to 914. If the classical objective function gradient with finite difference methods were used, each iteration of the optimization algorithm for the hetero/iso/linear-elastic forward elasticity problem would have to be solved 914 × 2 = 1828 times based on the central difference method to approximate the sensitivity of the objective function. As such, the optimization algorithm would be rendered prohibitively expensive and therefore, adjoint methods for optimization were utilized (Oberai et al., 2003). This effective method required solving the problem only twice; once for the forward elasticity problem and once for the adjoint problem.
The derivation of the adjoint method to formulate our IFEA problem employed the objective functional below, (with the predefined measure of error between simulation and experiment):where was the objective functional, vector was the global displacement vector, and was the set of unknown material parameters. Note that was dependent on the because knowing the material parameters allowed us to run the forward elasticity problem and solve for the nodal displacements. Therefore, the dependence of to was implicit. The size of the vectors u and p were and , where was the total degrees of the freedom of the FE model. Specifically, our mesh had 263 nodes where each node has three translational degrees of freedom (no rotational degrees), and therefore, . The size of vector for 457 meshed elements with two unknown material parameters and was . Since represented a measure of error, it took the following form:where the first term was the L-2 norm of the error, and the second term as a regularization parameter to tackle the ill-posed aspect of the inverse problem (Isakov, 2017). One popular choice for could be , where was a regularization parameter selected to be a very small number (10−6) selected based on the theory of residues (Tikhonov et al., 2013). The regularization was only applied to the heterogeneous model with the gradient-based optimization. In our case, referred to the five time points at which the error was calculated through the full inflation path.
The derivative of for the optimization algorithm was defined as:
The values for and were known given their explicit definition in Eq. 1. In order to get the term , the forward elasticity problem had to be solved since in general we do not have an analytical relation between and . The forward elasticity problem was cast into the following standard discretized format obtained from the FE model:where was the global stiffness matrix, and was the global load vector. From there the partial derivative of Eq. 4 was:whereor in the index notation
Substituting back into Eq. 3 yielded
The only problematic term was because it required solving the forward elasticity problem each time a small change was made to our unknown parameter . To address this issue, we wrotewhere was named the adjoint variable. Rearranging yielded
Note that from the definition, and from the symmetry of stiffness matrix. Therefore, to solve the adjoint equation, we applied the difference between the simulation and experiment displacements to drive the forward elasticity problem. Having solved for , the gradient was written aswhere acted as a Lagrange multiplier. For our simpler case with no regularization and the external load being independent of the unknown material properties , the derivative of the objective functional was simplified to
To calculate this simplified objective functional gradient, the FE problem was solved two times: one was the forward elasticity problem to obtain and the second one was to solve the adjoint set of equations to get . To obtain the term we slightly perturbed the material property and collected the assembled stiffness matrix by using the matrix generation procedure available in Abaqus.
Particle Swarm Optimization Algorithm
In this well-known algorithm (Kennedy and Eberhart, 1995) a random population of nPop particles was initially generated. Each particle was basically a point in our search space for the optimal material properties. For example, in the homo/iso/hyper case, the position of each particle was defined by its value of and randomly drawn from a specified range given in Table 1. The FE model for each particle was solved and the corresponding error was evaluated. In order to update the position of the particles toward the location of the global minima, the velocity of each particle was updated based on:
This determined the direction along which the value of the material properties was to be changed (i.e., increased or decreased). In Eq. 9, was a damping factor which reduced the momentum of the particles as they iteratively progressed towards finding the global optima (Kennedy and Eberhart, 1995); it started from 1.0 and was multiplied by a constant of 0.99 after each iteration. The parameters c1 and c2 controlled the local and global search weights, respectively. and referred to the best score (smallest error) of a given particle throughout the whole iterations passed thus far and the one within the current iteration, respectively. The parameter referred to the best score of the overall population. The particles position was then updated by adding the calculated velocity to the current position. FE simulations were then performed for the whole batch of particles in an iterative fashion until the optimization algorithm converged.
In implementing the PSO algorithm, it was important to impose proper upper and lower bounds on the velocity and position vectors of each particle to avoid local minima or particles getting stuck in the neighborhood of each other or worse yet, on the boundaries (Kennedy and Eberhart, 1995). Prior to updating the position vectors using Eq. 9, any element of the velocity vector that had values above or below were set to or , respectively. Then we checked for the position vectors: if a particle’s position was beyond the limits defined by the range , we checked its velocity vector; if it was pointing outside the position bound (meaning adding the velocity to the position would have resulted in the particle position landing outside of its permitted range), velocity component was set to zero, and the particle position was set to the corresponding upper or lower bound. The parameters and in Eq. 9 were set to be 2.0 (Kennedy and Eberhart, 1995). The ranges for particle velocity were set as for the homo/iso/hyper and homo/aniso/hyper cases, and as for the hetero/iso/linear-elastic case. The parameter nPop was set as 24, 48, and 1,000 for the homo/hyper/iso, homo/aniso/hyper, and hetero/iso/linear-elastic cases, respectively. These hyperparameters maintained a wide parameter value range, helped the algorithm converge better, and were tuned based on our preliminary sensitivity analysis studies.
Results
The Interpolated Deformation and Strain Measures of the Lung
The experimental displacement values were imposed on the FE model using the generic homo/iso/hyper case to confirm the validity of the interpolation technique and the strain orientations against the DIC system calculations. The data matched nearly identically for each of the five inflation stages (0.4, 0.8, 1.2, 1.6, and 2.0 s) as shown in Figure 3, and was subsequently used in the optimization scheme. The motion of the lung during inspiration was substantially inhomogeneous: the anterior region of the lobe exhibited the greatest distention and pronounced aeration (as much as 25 mm). The imposed displacements at the nodes along the lobe perimeter were also non-zero and were interpolated to represent the actual boundary conditions of the deforming surface properly. The maximum and minimum in-plane principal strains (major and minor strains) obtained at the full inflation stage were shown in Figure 3B and Figure 3C, valuing no more than 0.5 and 0.15, respectively. The major strain predominantly aligned with the medial-lateral direction while the minor strain was preferentially aligned with the anterior-posterior direction.
FIGURE 3
The Optimized Compound Material Properties of the Lung
The set of optimal material properties for the homo/iso/hyper and homo/aniso/hyper cases, and for a wide array of starting points, was listed in Tables 1, 2, respectively. The three material properties’ average ±standard deviation for the homo/iso/hyper was kPa, kPa, and kPa−1. The calculated shear and bulk moduli were kPa, and MPa, respectively. The parameter , which controlled the contribution of , consistently converged to its allowed lower bound of 1.0 kPa and was two order of magnitudes smaller than . The smallness of C01 in comparison to C10 indicates that the strain energy function is largely controlled by the stretch response of the tissue and not the distortion part. In ventilating the lung, the DIC strains also suggested that shear strains were minimal. Another consequence of C01 being very small is that the stress-strain curve would start to look less nonlinear.
The optimal set of material properties calculated for the homo/aniso/hyper case for seven optimization runs was kPa, kPa−1, kPa, , and (Table 2). The shear modulus was kPa, comparable to the 275 kPa obtained for the previous homo/iso/hyper case.
While the three material parameters , , and κ converged to their optimal values, the model did not depend on and . The anisotropic hyperelastic formulation in Abaqus ignores the compressibility coefficient and therefore, the objective function was simply insensitive to this parameter and D remained unchanged (Abaqus: Theory Manual, 2011). Conversely, k2, which was a dimensionless parameter exponentially controlling the contribution of the fibers in the overall strain energy function, did not yield a specific value because always converged to 0.33 and effectively zeroed out the anisotropic part of the strain energy density function. Therefore, , contributing to the anisotropic term, did not converge to a meaningful value because it played no role in the strain density function. Despite accounting for anisotropic lung behavior, the inverse optimization framework found no anisotropic advantage over the isotropic model. Given this observation, the homo/aniso/hyper case was not pursued further, and the considered models were limited to the homo/iso/hyper and hetero/iso/linear-elastic cases.
For the hetero/iso/linear-elastic case, we plotted the shear modulus map of the lung shown in Figure 4. In both SQP and PSO optimization schemes, the resulting spatial distribution of the shear modulus exhibited strong heterogeneity. The value for shear and bulk moduli was 108–312 kPa and 144 kPa–17.2 MPa, respectively, with the tissue softening from the posterior to anterior regions of the lobe.
FIGURE 4
The displacements and strains measured by DIC and predicted by the homogeneous and heterogeneous IFEA model were shown in Figure 5. The predictions of the heterogeneous model better matched the DIC displacement fields compared to that of the homogenous model. The overinflation of the lung at the anterior region was particularly well predicted by the heterogeneous model (Figure 5B).
FIGURE 5
The three components of displacement errors (percent normalized with respect of the maximum tissue displacement) for homogeneous and heterogeneous cases were shown in Figure 6. The errors for the heterogeneous case were consistently smaller than that of the homogeneous case and greatest in the z-direction (i.e., the direction at which DIC camera overlooked the lung), likely corresponding to the largest displacement values also being in the z-direction.
FIGURE 6
The overall IFEA settings and results of the two optimization algorithms applied to the homo/iso/hyper and hetero/iso/linear-elastic cases were summarized in Table 3. The homogeneous isotropic model returned an average error of 2.3 mm for both gradient-based and non-gradient-based optimization schemes. The average error for the heterogeneous lung model was consistently smaller than the homogeneous case; the PSO algorithm resulted in a slightly smaller error compared to the SQP algorithm (1.3 vs. 1.6 mm). The relative computational cost (CPU time normalized with respect to the fastest case) and number of iterations to reach the optimal solution were also given; the PSO algorithms for the homo/iso/hyper and hetero/iso/linear-elastic models were the fastest and most expensive simulations, respectively.
TABLE 3
| Homo/Iso/Hyper | Hetero/Iso/Linear-elastic | ||
| Gradient-based optimization | Error (mm) | 2.3 | 1.6 |
| Relative cost | 2.7 | 3.1 | |
| Algorithm | TRR | SQP | |
| Iterations | 17 | 7 | |
| Non-gradient-based optimization | Error (mm) | 2.3 | 1.3 |
| Relative cost | 1.0 | 68.0 | |
| Algorithm | PSO | PSO | |
| Iterations | 30 | 77 | |
General IFEA settings and results.
Discussion
While a one-to-one comparison between our reduced-order model and extracted pulmonary tissue specimen measures is impractical, we strived to compare the compound material parameters (averaging shear modulus of hundreds of kPa in all three constitutive models) with the reported ranges for parenchyma, airway, and pleura layer individual component responses (Tables 1, 2 and Figure 4). Our shear modulus values, representing the combined parenchyma, airway, and pleura layer materials, are greater than those of isolated lung parenchyma (0.17–0.27 kPa) and the alveolar wall (1.74 kPa), comparatively estimated by converting the previously reported elastic modulus with an assumed Poisson ratio of 0.43 (Lai-Fook et al., 1976; Cavalcante et al., 2005). Similarly converting the reported 10–70 kPa elastic moduli of isolated airway specimens (Eskandari et al., 2019) yields airway shear modulus range of 3.5–25 kPa, also less than our combined shear modulus results. Conversely, the encapsulating visceral pleura layer is approximated to have a shear modulus of ∼200 kPa at low stretch ratios (Humphrey et al., 1986), similar to our values.
The ratio between the bulk and shear modulus (often used to gage material compressibility) was nearly 5.5 for the homo/iso/hyper case, and 0.5 for the hetero/iso/linear-elastic case. This is significantly less than incompressible materials (with a ratio greater than 1,000) and justifies the use of the compressible material model, as previously suggested (Birzle and Wall, 2019). Our model suggests that lung elasticity is not distributed evenly across the regions, but that the shear moduli is smallest in the anterior region (Figure 4), corresponding to the location of maximum deformation. While literature substantiating this tissue heterogeneity across the organ is not yet available, regionally extracted tissue subjected to tensile or indentation tests can enable future comparisons.
Using this reduced-order model of the lung, our optimization scheme finds the anisotropic material model can be interchanged with an isotropic representation (since = 0.33) and still sufficiently capture the experimental displacements. While this finding is bound to the model limitation and calls for further experimental works to validate, it still can be substantiated given the isotropic material behavior of the parenchyma (Fung, 1988) where collagen and elastin fibers are randomly oriented (Toshima et al., 2004). However, the major strains were found to predominantly align with the medial-lateral direction while the minor strains were preferentially aligned with the anterior-posterior direction (Figure 3); this indicates that the spatial patterns and strain orientations were possibly a result of the geometry and loading of the lung more so than the anisotropic nature of the tissue material itself. An alternative hypothesis is the embedded monopodial main bronchial airway, which delivers oxygen from the anterior to posterior region and is twice as compliant circumferentially than axially, enables greater stretch in the medial-lateral direction (Sattari and Eskandari, 2020). Therefore, it is plausible that larger collagen-enriched airways may contribute to the anisotropic strain distribution in the lung to a great extent. Including a model mapping of the major airway pathways may help further differentiate the tissue matrix versus the effect of the structural reinforcement. This hypothesis will be explored further in upcoming mice and human lung experiments with differing bronchial branching patterns and collateral ventilation compared to the pig.
This reduced-order in-silico model of the lung facilitates a novel and much-needed class of inverse modeling approaches for the respiratory system. Current pulmonary biomechanical models of the lung can be categorized into two classes: 1) models primarily based on in-vivo kinematics data obtained from computed tomography (CT) or magnetic resonance (MR) images (Al-Mayah et al., 2010; Eom et al., 2010; Li et al., 2013; Ilegbusi et al., 2016; Ladjal et al., 2021), and 2) classic models idealizing the lung as single/multi resistive compartments calibrated with pressure-volume data (Bates, 2009). While these methods have been quite insightful, their shortcomings have motivated the novel approach put forward in this study. For instance, CT- and MR-based models (class 1) utilize convoluted and tedious deformable image registration (DIR), further challenged by the lack of a universal ground truth of lung nodules which necessitates expert-determined anatomical landmark detection and hinders model validation (Sotiras et al., 2013; Sarrut et al., 2017). Additionally, compartmentalized models (class 2) describe global bulk elastance and resistance behaviors and neglect the intricate multiscale architecture of the lung, omitting the local heterogeneities and strain risers responsible for inflammation and damage (Vlahakis et al., 1999; Gattinoni et al., 2003). The absence of controllable testing parameters and continuous measures in-vivo limits basic lung kinetic and kinematic investigations, such as exploring the role of ventilation volume and rate, and contrasting physiological negative-versus artificial positive-pressure ventilation (Eskandari et al., 2021). Our ex-vivo informed in-silico approach can readily establish the mechanical science of breathing by merging detailed data acquisition with computational predictions.
Forging a bridge between tissue-scale kinematics and organ-level kinetics facilitates fundamental explorations of multiscale characterization and inaugurates several applications. Generalizing and informing the model with multiple lung data sets and complete inhalation and exhalation breathing pressures can empower surgical planning strategies based on minimizing changes in strain patterns from lobectomies, segmentectomies, and wedge resections to preoperatively improve patient outcomes instead of postoperative evaluations (Charloux and Quoix, 2017). Extending this framework to include diseased lungs categories, such as fibrosis or asthma, will enable regional tissue remodeling detection studies to discover load- and deformation-based pathological deviations to guide therapies, as inspired by similar FE-based models of the lung being used to better aerosol deposition in asthmatic patients (Wall et al., 2010; Soni and Aliabadi, 2013) and optimal radiation therapy (Werner et al., 2009; Ilegbusi et al., 2016). Furthermore, this model has potential applications for the design of ventilators: many didactic (Anderson et al., 2009) or clinical ventilators (Zuckerberg et al., 2020) contain a lung-replicating rubber bladder or elastic balloon component with mechanical properties that are not physiologically representative. Our FE lung membrane-like model and optimized material properties can enhance the design resemblance of these ventilator systems to improve patient care.
This study has several limitations. Firstly, one animal was used to demonstrate this IFEA framework and therefore, statistically conclusive results regarding the material property values are not warranted. Second, while the optimization algorithms minimized the error, the final error was still not completely vanished; this indicates lung-specific constitutive models are needed to better replicate the DIC measurements, as concluded by earlier works (Eskandari et al., 2019). Third, while DIC allows continuous and evolutionary behaviors of the lung to be examined whereas digital volume correlation techniques are at discrete snapshots (Arora et al., 2021), DIC can only access the lung surface and the internal structure of the lung and the volumetric strain distributions are not represented; a potential enhancement to this technique could be the use of mirrors and prisms to collect multi-angled views. Fourth, the framework is built on ex-vivo setting, hence the predicted material properties are likely to be influenced by the deformation and kinetics of the ribcage and diaphragm. Fifth, the anisotropic formulation implemented in the FE-package used in this study does not take into account the compressibility aspect of plane-stress elements, such as membrane elements in our model and assumes constant volume instead while still allowing for reduction in the element thickness due to in-plane deformations. Thus, the compressibility coefficient (parameter D in the HGO model) is not optimized by the IFEA framework. New implementation of the HGO formulation that do indeed account for the compressibility are needed to further investigate the effects of this limitation and provide insight into the compressibility of the anisotropic model. Lastly, for simplicity and computational time considerations, the heterogeneous model was a linear model whereas a nonlinear model might result in a better calibrated model.
Conclusion
This study establishes a computational model representing local lung kinetics by associating global organ-level pressures and volumes to tissue-level kinematics. This is achieved by developing a novel lung application IFEA framework informed and verified by ex-vivo continuous DIC measurements from a porcine lung controlled via a custom-designed respiration apparatus. The resulting FE model introduces a model constructed solely from the geometry and deformation of its external surface as a result of the applied inflation load. This in-silico reduced-order pulmonary surrogate consolidates complex lung tissues (i.e., the visceral pleura, bulk parenchymal tissue, and the airway tree) into a simplified 3-D surface model, yielding compound material properties of a membrane representative of pressure-deformation features of the lung. Furthermore, the heterogeneous lung elasticity map presented in this study empowers new avenues to improve characterization of diseased states by enabling region-specific assessments of mechanical remodeling, such as variations in tissue elasticity, thus far critically absent in the field.
Statements
Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.
Ethics statement
Ethical review and approcal was not required for the specimens used in the study because the animal lungs were obtained from an abattoir.
Author contributions
MM-G designed the model and optimization schemes. CM and SS performed experiments. ME conceptualized and supervised research. ME and MM-G interpreted results and analysed the data. HA, ME, and MM-G wrote the manuscript. All authors contributed to the article and approved the submitted version.
Funding
The authors gratefully acknowledge funding support from the Dassault Systèmes U.S. Foundation Grant and the Hellman Fellows Program to ME, the Sêr Cymru programme, Welsh Government for supporting HA, and the Global Wales International Research Mobility Fund (UNIW/RMF-SU/03) for facilitating the transatlantic collaboration. This work is dedicated to the millions who lost their lives due to COVID-19.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1
Abaqus Theory Manual (2011). Abaqus: Theory Manual. Providence, RI, USA: Dassault Systèmes Simulia Corp.
2
Al-MayahA.MoseleyJ.VelecM.HunterS.BrockK. (2010). Deformable Image Registration of Heterogeneous Human Lung Incorporating the Bronchial Tree. Med. Phys.37, 4560–4571. 10.1118/1.3471020
3
AltmanN. S. (1992). An Introduction to Kernel and Nearest-Neighbor Nonparametric Regression. The Am. Statistician46, 175–185. 10.1080/00031305.1992.10475879
4
AmatoM. B. P.MeadeM. O.SlutskyA. S.BrochardL.CostaE. L. V.SchoenfeldD. A.et al (2015). Driving Pressure and Survival in the Acute Respiratory Distress Syndrome. N. Engl. J. Med.372, 747–755. 10.1056/NEJMsa1410639
5
AndersonJ.GoplenC.MurrayL.SeashoreK.SoundarrajanM.LokutaA.et al (2009). Human Respiratory Mechanics Demonstration Model. Adv. Physiol. Educ.33, 53–59. 10.1152/advan.90177.2008
6
AroraH.NilaA.VitharanaK.SherwoodJ. M.NguyenT.-T. N.KarunaratneA.et al (2017). Microstructural Consequences of Blast Lung Injury Characterized with Digital Volume Correlation. Front. Mater.4, 41. 10.3389/fmats.2017.00041
7
AroraH.MitchellR.JohnstonR.ManolesosM.HowellsD.SherwoodJ.et al (2021). Correlating Local Volumetric Tissue Strains with Global Lung Mechanics Measurements. Materials14, 439. 10.3390/ma14020439
8
AtkesonA. (2020). What Will Be the Economic Impact of COVID-19 in the US? Rough Estimates of Disease Scenarios, (No. w26867). Natl. Bur. Econ. Res.10.3386/w26867
9
BaiT. R.KnightD. A. (2005). Structural Changes in the Airways in Asthma: Observations and Consequences. Clin. Sci.108, 463–477. 10.1042/CS20040342
10
BatesJ. H. T. (2009). Lung Mechanics: An Inverse Modeling Approach. Cambridge University Press.
11
BirzleA. M.WallW. A. (2019). A Viscoelastic Nonlinear Compressible Material Model of Lung Parenchyma - Experiments and Numerical Identification. J. Mech. Behav. Biomed. Mater.94, 164–175. 10.1016/j.jmbbm.2019.02.024
12
BirzleA. M.HobrackS. M. K.MartinC.UhligS.WallW. A. (2019). Constituent-Specific Material Behavior of Soft Biological Tissue: Experimental Quantification and Numerical Identification for Lung Parenchyma. Biomech. Model. Mechanobiol.18, 1383–1400. 10.1007/s10237-019-01151-3
13
BoggsP. T.TolleJ. W. (1995). Sequential Quadratic Programming. Acta Numerica4, 1–51. 10.1017/S0962492900002518
14
BoyceB. L.GrazierJ. M.JonesR. E.NguyenT. D. (2008). Full-Field Deformation of Bovine Cornea under Constrained Inflation Conditions. Biomaterials29, 3896–3904. 10.1016/j.biomaterials.2008.06.011
15
CavalcanteF. S. A.ItoS.BrewerK.SakaiH.AlencarA. M.AlmeidaM. P.et al (2005). Mechanical Interactions between Collagen and Proteoglycans: Implications for the Stability of Lung Tissue. J. Appl. Physiol.98, 672–679. 10.1152/japplphysiol.00619.2004
16
Centers for Disease Control and Prevention (2015). Health U. S. 2015 US Dep. Health Hum. Serv. Cent. Dis. Control. Prev. Natl. Cent. Health Stat. DHHS Publ., 2016–1232.
17
CharlouxA.QuoixE. (2017). Lung Segmentectomy: Does it Offer a Real Functional Benefit over Lobectomy?Eur. Respir. Rev.26, 170079. 10.1183/16000617.0079-2017
18
ChuT. C.RansonW. F.SuttonM. A. (1985). Applications of Digital-Image-Correlation Techniques to Experimental Mechanics. Exp. Mech.25, 232–244. 10.1007/BF02325092
19
CignoniP.CallieriM.CorsiniM.DellepianeM.GanovelliF.RanzugliaG. (2008). “MeshLab: An Open-Source Mesh Processing Tool,” in Eurographics Ital. Chapter Conf. 2008, 129–136.
20
DreyfussD.SaumonG. (1998). Ventilator-Induced Lung Injury. Am. J. Respir. Crit. Care Med.157, 294–323. 10.1164/ajrccm.157.1.9604014
21
EomJ.XuX. G.DeS.ShiC. (2010). Predictive Modeling of Lung Motion Over the Entire Respiratory Cycle Using Measured Pressure-Volume Data, 4DCT Images, and Finite-Element Analysis. Med. Phys.37, 4389–4400. 10.1118/1.3455276
22
EskandariM.KuhlE. (2015). Systems Biology and Mechanics of Growth. Wires Syst. Biol. Med.7, 401–412. 10.1002/wsbm.1312
23
EskandariM.PfallerM.KuhlE. (2013). On the Role of Mechanics in Chronic Lung Disease. Materials6, 5639–5658. 10.3390/ma6125639
24
EskandariM.KuschnerW. G.KuhlE. (2015). Patient-Specific Airway Wall Remodeling in Chronic Lung Disease. Ann. Biomed. Eng.43, 2538–2551. 10.1007/s10439-015-1306-7
25
EskandariM.JaviliA.KuhlE. (2016). Elastosis during Airway wall Remodeling Explains Multiple Co-Existing Instability Patterns. J. Theor. Biol.403, 209–218. 10.1016/j.jtbi.2016.05.022
26
EskandariM.ArvayoA. L.LevenstonM. E. (2018). Mechanical Properties of the Airway Tree: Heterogeneous and Anisotropic Pseudoelastic and Viscoelastic Tissue Responses. J. Appl. Physiol.125, 878–888. 10.1152/japplphysiol.00090.2018
27
EskandariM.NordgrenT. M.O’ConnellG. D. (2019). Mechanics of Pulmonary Airways: Linking Structure to Function through Constitutive Modeling, Biochemistry, and Histology. Acta Biomater.97, 513–523. 10.1016/j.actbio.2019.07.020
28
EskandariM.SattariS.MarianoC. A. (2021). “Investigating the Mechanics of Positive- versus Negative-Pressure Ventilation,” in TP127 LUNG INJURY AND MECHANICAL VENTILATION, American Thoracic Society International Conference Abstracts (American Thoracic Society), A4671. 10.1164/ajrccm-conference.2021.203.1_MeetingAbstracts.A4671
29
FaffeD. S.ZinW. A. (2009). Lung Parenchymal Mechanics in Health and Disease. Physiol. Rev.89, 759–775. 10.1152/physrev.00019.2007
30
FungY. C. (1988). A Model of the Lung Structure and its Validation. J. Appl. Physiol.64, 2132–2141. 10.1152/jappl.1988.64.5.2132
31
GattinoniL.CarlessoE.CadringherP.ValenzaF.VagginelliF.ChiumelloD. (2003). Physical and Biological Triggers of Ventilator-Induced Lung Injury and its Prevention. Eur. Respir. J.22, 15s–25s. 10.1183/09031936.03.00021303
32
HolzapfelG. A.GasserT. C.OgdenR. W. (2000). A New Constitutive Framework for Arterial Wall Mechanics and a Comparative Study of Material Models. J. Elast. Phys. Sci. Sol.61, 1–48. 10.1023/A:1010835316564
33
HumphreyJ. D.VawterD. L.VitoR. P. (1986). Mechanical Behavior of Excised Canine Visceral Pleura. Ann. Biomed. Eng.14, 451–466. 10.1007/BF02367365
34
SeyfiB.SanthanamA. P.IlegbusiO. J. (2016). A Biomechanical Model of Human Lung Deformation Utilizing Patient-Specific Elastic Property. J. Cancer Ther.07, 402–415. 10.4236/jct.2016.76043
35
IsakovV. (2017). “Inverse Problems for Partial Differential Equations,” in Applied Mathematical Sciences. 3rd ed (Springer). 10.1007/978-3-319-51658-5
36
JonesE. M.IadicolaM. A. (2018). A Good Practices Guide for Digital Image Correlation. Int. Digit. Image Correl. Soc.10. 10.32720/idics/gpg.ed1
37
KennedyJ.EberhartR. (1995). “Particle Swarm Optimization,” in Proceedings of the IEEE International Conference on Neural Networks, 1942–1948. 10.1109/ICNN.1995.488968
38
LadjalH.BeuveM.GiraudP.ShariatB. (2021). Towards Non-Invasive Lung Tumor Tracking Based on Patient Specific Model of Respiratory System. IEEE Trans. Biomed. Eng.68, 2730–2740. 10.1109/TBME.2021.3053321
39
Lai-FookS. J.WilsonT. A.HyattR. E.RodarteJ. R. (1976). Elastic Constants of Inflated Lobes of Dog Lungs. J. Appl. Physiol.40, 508–513. 10.1152/jappl.1976.40.4.508
40
LiM.CastilloE.ZhengX.-L.LuoH.-Y.CastilloR.WuY.et al (2013). Modeling Lung Deformation: A Combined Deformable Image Registration Method with Spatially Varying Young's Modulus Estimates. Med. Phys.40, 081902. 10.1118/1.4812419
41
Maghsoudi-GanjehM.SattariS.EskandariM. (2021). Mechanical Behavior of the Airway wall in Respiratory Disease. Curr. Opin. Physiol.22, 100445. 10.1016/j.cophys.2021.05.008
42
MallettK. F.ArrudaE. M. (2017). Digital Image Correlation-Aided Mechanical Characterization of the Anteromedial and Posterolateral Bundles of the Anterior Cruciate Ligament. Acta Biomater.56, 44–57. 10.1016/j.actbio.2017.03.045
43
MarianoC. A.SattariS.Maghsoudi-GanjehM.TartibiM.LoD. D.EskandariM. (2020). Novel Mechanical Strain Characterization of Ventilated Ex Vivo Porcine and Murine Lung Using Digital Image Correlation. Front. Physiol.11, 1536. 10.3389/fphys.2020.600492
44
MooneyM. (1940). A Theory of Large Elastic Deformation. J. Appl. Phys.11, 582–592. 10.1063/1.1712836
45
NocedalJ.StephenW. (2006). Numerical Optimization. Springer Science & Business Media.
46
OberaiA. A.GokhaleN. H.Feij oG. R. (2003). Solution of Inverse Problems in Elasticity Imaging Using the Adjoint Method. Inverse Probl.19, 297–313. 10.1088/0266-5611/19/2/304
47
SarrutD.BaudierT.AyadiM.TanguyR.RitS. (2017). Deformable Image Registration Applied to Lung SBRT: Usefulness and Limitations. Physica Med.44, 108–112. 10.1016/j.ejmp.2017.09.121
48
SattariS.EskandariM. (2020). Characterizing the Viscoelasticity of Extra- and Intra-Parenchymal Lung Bronchi. J. Mech. Behav. Biomed. Mater.110, 103824. 10.1016/j.jmbbm.2020.103824
49
SattariS.MarianoC. A.VittalbabuS.VelazquezJ. V.PostmaJ.HorstC.et al (2020). Introducing a Custom-Designed Volume-Pressure Machine for Novel Measurements of Whole Lung Organ Viscoelasticity and Direct Comparisons between Positive- and Negative-Pressure Ventilation. Front. Bioeng. Biotechnol.8, 1183. 10.3389/fbioe.2020.578762
50
SoniB.AliabadiS. (2013). Large-Scale CFD Simulations of Airflow and Particle Deposition in Lung Airway. Comput. Fluids88, 804–812. 10.1016/j.compfluid.2013.06.015
51
SotirasA.DavatzikosC.ParagiosN. (2013). Deformable Medical Image Registration: A Survey. IEEE Trans. Med. Imaging32, 1153–1190. 10.1109/TMI.2013.2265603
52
SteihaugT. (1983). The Conjugate Gradient Method and Trust Regions in Large Scale Optimization. SIAM J. Numer. Anal.20, 626–637. 10.1137/0720042
53
SukiB.BatesJ. H. T. (2011). Lung Tissue Mechanics as an Emergent Phenomenon. J. Appl. Physiol.110, 1111–1118. 10.1152/japplphysiol.01244.2010
54
SukiB.LutchenK. R.IngenitoE. P. (2003). On the Progressive Nature of Emphysema. Am. J. Respir. Crit. Care Med.168, 516–521. 10.1164/rccm.200208-908PP
55
SuttonM. A.KeX.LessnerS. M.GoldbachM.YostM.ZhaoF.et al (2008). Strain Field Measurements on Mouse Carotid Arteries Using Microscopic Three-Dimensional Digital Image Correlation. J. Biomed. Mater. Res.84A, 178–190. 10.1002/jbm.a.31268
56
TawhaiM. H.BatesJ. H. T. (2011). Multi-scale Lung Modeling. J. Appl. Physiol.110, 1466–1472. 10.1152/japplphysiol.01289.2010
57
The Acute Respiratory Distress Syndrome Network (2000). Ventilation with Lower Tidal Volumes as Compared with Traditional Tidal Volumes for Acute Lung Injury and the Acute Respiratory Distress Syndrome. N. Engl. J. Med.342 (18), 1301–1308. 10.1056/NEJM200005043421801
58
TikhonovA. N.GoncharskyA.StepanovV. V.YagolaA. G. (2013). Numerical Methods for the Solution of Ill-Posed Problems. Springer Science & Business Media.
59
ToshimaM.OhtaniY.OhtaniO. (2004). Three-Dimensional Architecture of Elastin and Collagen Fiber Networks in the Human and Rat Lung. Arch. Histology Cytol.67, 31–40. 10.1679/aohc.67.31
60
VlahakisN. E.SchroederM. A.LimperA. H.HubmayrR. D. (1999). Stretch Induces Cytokine Release by Alveolar Epithelial Cells In Vitro. Am. J. Physiology-Lung Cell Mol. Physiol.277, L167–L173. 10.1152/ajplung.1999.277.1.L167
61
WallW. A.WiechertL.ComerfordA.RauschS. (2010). Towards a Comprehensive Computational Model for the Respiratory System. Int. J. Numer. Meth. Biomed. Engng.26, 807–827. 10.1002/cnm.1378
62
WernerR.EhrhardtJ.SchmidtR.HandelsH. (2009). Patient-Specific Finite Element Modeling of Respiratory Lung Motion Using 4D CT Image Data. Med. Phys.36, 1500–1511. 10.1118/1.3101820
63
ZengY. J.YagerD.FungY. C. (1987). Measurement of the Mechanical Properties of the Human Lung Tissue. J. Biomech. Eng.109, 169–174. 10.1115/1.3138661
64
ZuckerbergJ.ShaikM.WidmeierK.KilbaughT.NelinT. D. (2020). A Lung for All: Novel Mechanical Ventilator for Emergency and Low-Resource Settings. Life Sci.257, 118113. 10.1016/j.lfs.2020.118113
Summary
Keywords
lung mechanics, inverse finite element analysis, digital image correlation (DIC), heterogeneity, anisotropy, hyperelasticity, in-silico ventilation
Citation
Maghsoudi-Ganjeh M, Mariano CA, Sattari S, Arora H and Eskandari M (2021) Developing a Lung Model in the Age of COVID-19: A Digital Image Correlation and Inverse Finite Element Analysis Framework. Front. Bioeng. Biotechnol. 9:684778. doi: 10.3389/fbioe.2021.684778
Received
24 March 2021
Accepted
04 October 2021
Published
26 October 2021
Volume
9 - 2021
Edited by
Jerome Noailly, Pompeu Fabra University, Spain
Reviewed by
Abdelwahed Barkaoui, International University of Rabat, Morocco
Behzad Shariat, Université Claude Bernard Lyon 1, France
Thiranja Prasad Babarenda Gamage, The University of Auckland, New Zealand
Updates
Copyright
© 2021 Maghsoudi-Ganjeh, Mariano, Sattari, Arora and Eskandari.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Mona Eskandari, mona.eskandari@ucr.edu
This article was submitted to Biomechanics, a section of the journal Frontiers in Bioengineering and Biotechnology
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.