# Developing a Lung Model in the Age of COVID-19: A Digital Image Correlation and Inverse Finite Element Analysis Framework

^{1}Department of Mechanical Engineering, University of California, Riverside, Riverside, CA, United States^{2}Faculty of Science and Engineering, Swansea University, Swansea, United Kingdom^{3}BREATHE Center, School of Medicine, University of California, Riverside, Riverside, CA, United States^{4}Department of Bioengineering, University of California, Riverside, Riverside, CA, United States

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**. **(A)** The right lung lobe was selected for analysis and DIC strains are shown. **(B)** The pressure-volume data of the stabilized inflation cycle was used to extract the **(C)** pressure-time data applied to the FE model.

### 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**. The overall algorithm for the IFEA framework implementation. The model started with an initial estimate for the material parameters, then the optimization algorithm incrementally moved forward toward the optimal solution by minimizing the displacement error between model predictions and DIC measurements. In this workflow, u refers to the nodal displacement, σ is the Cauchy stress, ε is the technical strain.

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

##### Homogenous Anisotropic Hyperelastic Case

The Holzapfel-Gasser-Ogden (HGO) formulation (Holzapfel et al., 2000) with the strain energy density defined as *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, *x*-direction for each element. It should be noted that even though we define *x*-axis for the main fiber directions, since the parameter

##### Heterogeneous Isotropic Linear-Elastic Case

In this case we considered a linear elastic isotropic material model (Eskandari and Kuhl, 2015) with Young’s modulus

### 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**. 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**. 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 **u** and **p** were

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 ^{−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,

The derivative of

The values for

where

where

or in the index notation

Substituting

The only problematic term was

where

Note that

where

To calculate this simplified objective functional gradient, the FE problem was solved two times: one was the forward elasticity problem to obtain

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

This determined the direction along which the value of the material properties was to be changed (i.e., increased or decreased). In Eq. 9, *c*_{1} and *c*_{2} controlled the local and global search weights, respectively.

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 *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**. DIC measured displacements data, extracted at ∼7,000 evenly distributed points across the parenchymal surface were interpolated and applied to the FE model. **(A)** Displacement magnitude contour maps corresponding to five increments of 20% inflation steps. Anterior regions of the lung exhibited the largest distensions. Vector field map of major **(B)** and minor **(C)** strains were obtained from imposing interpolated DIC displacements to the model.

### 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 ^{−1}. The calculated shear and bulk moduli were

The optimal set of material properties calculated for the homo/aniso/hyper case for seven optimization runs was ^{−1},

While the three material parameters *D* remained unchanged (Abaqus: Theory Manual, 2011). Conversely, *k*_{2}, 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

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 shear modulus map obtained for the hetero/iso/linear-elastic case obtained from **(A)** the gradient-based and **(B)** derivative-free optimization algorithms. Both methods consistently demonstrated strong heterogeneity in the tissue elasticity distribution across the lung surface.

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**. Comparison between the displacement **(A–C)** and major strain **(D–F)** of the DIC and IFEA with either homogeneous and heterogeneous material models. The computed displacement and strain contours of the heterogeneous model agreed with the DIC data better than the homogeneous model. The heterogeneous model results are based on the SQP algorithm, and the PSO algorithm yielded similar results.

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 displacement error (Cartesian components) of the IFEA predictions assuming **(A)** heterogeneous and **(B)** homogeneous material models. The errors for the heterogenous model were consistently smaller. SQP algorithm shown for the heterogeneous model error and the PSO algorithm yielded similar results.

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.

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

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.

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

Abaqus Theory Manual (2011). *Abaqus: Theory Manual*. Providence, RI, USA: Dassault Systèmes Simulia Corp.

Al-Mayah, A., Moseley, J., Velec, M., Hunter, S., and Brock, K. (2010). Deformable Image Registration of Heterogeneous Human Lung Incorporating the Bronchial Tree. *Med. Phys.* 37, 4560–4571. doi:10.1118/1.3471020

Altman, N. S. (1992). An Introduction to Kernel and Nearest-Neighbor Nonparametric Regression. *The Am. Statistician* 46, 175–185. doi:10.1080/00031305.1992.10475879

Amato, M. B. P., Meade, M. O., Slutsky, A. S., Brochard, L., Costa, E. L. V., Schoenfeld, D. A., et al. (2015). Driving Pressure and Survival in the Acute Respiratory Distress Syndrome. *N. Engl. J. Med.* 372, 747–755. doi:10.1056/NEJMsa1410639

Anderson, J., Goplen, C., Murray, L., Seashore, K., Soundarrajan, M., Lokuta, A., et al. (2009). Human Respiratory Mechanics Demonstration Model. *Adv. Physiol. Educ.* 33, 53–59. doi:10.1152/advan.90177.2008

Arora, H., Nila, A., Vitharana, K., Sherwood, J. M., Nguyen, T.-T. N., Karunaratne, A., et al. (2017). Microstructural Consequences of Blast Lung Injury Characterized with Digital Volume Correlation. *Front. Mater.* 4, 41. doi:10.3389/fmats.2017.00041

Arora, H., Mitchell, R., Johnston, R., Manolesos, M., Howells, D., Sherwood, J., et al. (2021). Correlating Local Volumetric Tissue Strains with Global Lung Mechanics Measurements. *Materials* 14, 439. doi:10.3390/ma14020439

Atkeson, A. (2020). What Will Be the Economic Impact of COVID-19 in the US? Rough Estimates of Disease Scenarios, (No. w26867). *Natl. Bur. Econ. Res.* doi:10.3386/w26867

Bai, T. R., and Knight, D. A. (2005). Structural Changes in the Airways in Asthma: Observations and Consequences. *Clin. Sci.* 108, 463–477. doi:10.1042/CS20040342

Birzle, A. M., and Wall, W. A. (2019). A Viscoelastic Nonlinear Compressible Material Model of Lung Parenchyma - Experiments and Numerical Identification. *J. Mech. Behav. Biomed. Mater.* 94, 164–175. doi:10.1016/j.jmbbm.2019.02.024

Birzle, A. M., Hobrack, S. M. K., Martin, C., Uhlig, S., and Wall, W. A. (2019). Constituent-Specific Material Behavior of Soft Biological Tissue: Experimental Quantification and Numerical Identification for Lung Parenchyma. *Biomech. Model. Mechanobiol.* 18, 1383–1400. doi:10.1007/s10237-019-01151-3

Boggs, P. T., and Tolle, J. W. (1995). Sequential Quadratic Programming. *Acta Numerica* 4, 1–51. doi:10.1017/S0962492900002518

Boyce, B. L., Grazier, J. M., Jones, R. E., and Nguyen, T. D. (2008). Full-Field Deformation of Bovine Cornea under Constrained Inflation Conditions. *Biomaterials* 29, 3896–3904. doi:10.1016/j.biomaterials.2008.06.011

Cavalcante, F. S. A., Ito, S., Brewer, K., Sakai, H., Alencar, A. M., Almeida, M. P., et al. (2005). Mechanical Interactions between Collagen and Proteoglycans: Implications for the Stability of Lung Tissue. *J. Appl. Physiol.* 98, 672–679. doi:10.1152/japplphysiol.00619.2004

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.

Charloux, A., and Quoix, E. (2017). Lung Segmentectomy: Does it Offer a Real Functional Benefit over Lobectomy? *Eur. Respir. Rev.* 26, 170079. doi:10.1183/16000617.0079-2017

Chu, T. C., Ranson, W. F., and Sutton, M. A. (1985). Applications of Digital-Image-Correlation Techniques to Experimental Mechanics. *Exp. Mech.* 25, 232–244. doi:10.1007/BF02325092

Cignoni, P., Callieri, M., Corsini, M., Dellepiane, M., Ganovelli, F., and Ranzuglia, G. (2008). “MeshLab: An Open-Source Mesh Processing Tool,” in Eurographics Ital. Chapter Conf. 2008, 129–136.

Dreyfuss, D., and Saumon, G. (1998). Ventilator-Induced Lung Injury. *Am. J. Respir. Crit. Care Med.* 157, 294–323. doi:10.1164/ajrccm.157.1.9604014

Eom, J., Xu, X. G., De, S., and Shi, C. (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. doi:10.1118/1.3455276

Eskandari, M., and Kuhl, E. (2015). Systems Biology and Mechanics of Growth. *Wires Syst. Biol. Med.* 7, 401–412. doi:10.1002/wsbm.1312

Eskandari, M., Pfaller, M., and Kuhl, E. (2013). On the Role of Mechanics in Chronic Lung Disease. *Materials* 6, 5639–5658. doi:10.3390/ma6125639

Eskandari, M., Kuschner, W. G., and Kuhl, E. (2015). Patient-Specific Airway Wall Remodeling in Chronic Lung Disease. *Ann. Biomed. Eng.* 43, 2538–2551. doi:10.1007/s10439-015-1306-7

Eskandari, M., Javili, A., and Kuhl, E. (2016). Elastosis during Airway wall Remodeling Explains Multiple Co-Existing Instability Patterns. *J. Theor. Biol.* 403, 209–218. doi:10.1016/j.jtbi.2016.05.022

Eskandari, M., Arvayo, A. L., and Levenston, M. E. (2018). Mechanical Properties of the Airway Tree: Heterogeneous and Anisotropic Pseudoelastic and Viscoelastic Tissue Responses. *J. Appl. Physiol.* 125, 878–888. doi:10.1152/japplphysiol.00090.2018

Eskandari, M., Nordgren, T. M., and O’Connell, G. D. (2019). Mechanics of Pulmonary Airways: Linking Structure to Function through Constitutive Modeling, Biochemistry, and Histology. *Acta Biomater.* 97, 513–523. doi:10.1016/j.actbio.2019.07.020

Eskandari, M., Sattari, S., and Mariano, C. 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. doi:10.1164/ajrccm-conference.2021.203.1_MeetingAbstracts.A4671

Faffe, D. S., and Zin, W. A. (2009). Lung Parenchymal Mechanics in Health and Disease. *Physiol. Rev.* 89, 759–775. doi:10.1152/physrev.00019.2007

Fung, Y. C. (1988). A Model of the Lung Structure and its Validation. *J. Appl. Physiol.* 64, 2132–2141. doi:10.1152/jappl.1988.64.5.2132

Gattinoni, L., Carlesso, E., Cadringher, P., Valenza, F., Vagginelli, F., and Chiumello, D. (2003). Physical and Biological Triggers of Ventilator-Induced Lung Injury and its Prevention. *Eur. Respir. J.* 22, 15s–25s. doi:10.1183/09031936.03.00021303

Holzapfel, G. A., Gasser, T. C., and Ogden, R. 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. doi:10.1023/A:1010835316564

Humphrey, J. D., Vawter, D. L., and Vito, R. P. (1986). Mechanical Behavior of Excised Canine Visceral Pleura. *Ann. Biomed. Eng.* 14, 451–466. doi:10.1007/BF02367365

Seyfi, B., Santhanam, A. P., and Ilegbusi, O. J. (2016). A Biomechanical Model of Human Lung Deformation Utilizing Patient-Specific Elastic Property. *J. Cancer Ther.* 07, 402–415. doi:10.4236/jct.2016.76043

Isakov, V. (2017). “Inverse Problems for Partial Differential Equations,” in *Applied Mathematical Sciences*. 3rd ed (Springer). doi:10.1007/978-3-319-51658-5

Jones, E. M., and Iadicola, M. A. (2018). A Good Practices Guide for Digital Image Correlation. *Int. Digit. Image Correl. Soc.* 10. doi:10.32720/idics/gpg.ed1

Kennedy, J., and Eberhart, R. (1995). “Particle Swarm Optimization,” in Proceedings of the IEEE International Conference on Neural Networks, 1942–1948. doi:10.1109/ICNN.1995.488968

Ladjal, H., Beuve, M., Giraud, P., and Shariat, B. (2021). Towards Non-Invasive Lung Tumor Tracking Based on Patient Specific Model of Respiratory System. *IEEE Trans. Biomed. Eng.* 68, 2730–2740. doi:10.1109/TBME.2021.3053321

Lai-Fook, S. J., Wilson, T. A., Hyatt, R. E., and Rodarte, J. R. (1976). Elastic Constants of Inflated Lobes of Dog Lungs. *J. Appl. Physiol.* 40, 508–513. doi:10.1152/jappl.1976.40.4.508

Li, M., Castillo, E., Zheng, X.-L., Luo, H.-Y., Castillo, R., Wu, Y., et al. (2013). Modeling Lung Deformation: A Combined Deformable Image Registration Method with Spatially Varying Young's Modulus Estimates. *Med. Phys.* 40, 081902. doi:10.1118/1.4812419

Maghsoudi-Ganjeh, M., Sattari, S., and Eskandari, M. (2021). Mechanical Behavior of the Airway wall in Respiratory Disease. *Curr. Opin. Physiol.* 22, 100445. doi:10.1016/j.cophys.2021.05.008

Mallett, K. F., and Arruda, E. M. (2017). Digital Image Correlation-Aided Mechanical Characterization of the Anteromedial and Posterolateral Bundles of the Anterior Cruciate Ligament. *Acta Biomater.* 56, 44–57. doi:10.1016/j.actbio.2017.03.045

Mariano, C. A., Sattari, S., Maghsoudi-Ganjeh, M., Tartibi, M., Lo, D. D., and Eskandari, M. (2020). Novel Mechanical Strain Characterization of Ventilated *Ex Vivo* Porcine and Murine Lung Using Digital Image Correlation. *Front. Physiol.* 11, 1536. doi:10.3389/fphys.2020.600492

Mooney, M. (1940). A Theory of Large Elastic Deformation. *J. Appl. Phys.* 11, 582–592. doi:10.1063/1.1712836

Oberai, A. A., Gokhale, N. H., and Feij o, G. R. (2003). Solution of Inverse Problems in Elasticity Imaging Using the Adjoint Method. *Inverse Probl.* 19, 297–313. doi:10.1088/0266-5611/19/2/304

Sarrut, D., Baudier, T., Ayadi, M., Tanguy, R., and Rit, S. (2017). Deformable Image Registration Applied to Lung SBRT: Usefulness and Limitations. *Physica Med.* 44, 108–112. doi:10.1016/j.ejmp.2017.09.121

Sattari, S., and Eskandari, M. (2020). Characterizing the Viscoelasticity of Extra- and Intra-Parenchymal Lung Bronchi. *J. Mech. Behav. Biomed. Mater.* 110, 103824. doi:10.1016/j.jmbbm.2020.103824

Sattari, S., Mariano, C. A., Vittalbabu, S., Velazquez, J. V., Postma, J., Horst, C., 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. doi:10.3389/fbioe.2020.578762

Soni, B., and Aliabadi, S. (2013). Large-Scale CFD Simulations of Airflow and Particle Deposition in Lung Airway. *Comput. Fluids* 88, 804–812. doi:10.1016/j.compfluid.2013.06.015

Sotiras, A., Davatzikos, C., and Paragios, N. (2013). Deformable Medical Image Registration: A Survey. *IEEE Trans. Med. Imaging* 32, 1153–1190. doi:10.1109/TMI.2013.2265603

Steihaug, T. (1983). The Conjugate Gradient Method and Trust Regions in Large Scale Optimization. *SIAM J. Numer. Anal.* 20, 626–637. doi:10.1137/0720042

Suki, B., and Bates, J. H. T. (2011). Lung Tissue Mechanics as an Emergent Phenomenon. *J. Appl. Physiol.* 110, 1111–1118. doi:10.1152/japplphysiol.01244.2010

Suki, B., Lutchen, K. R., and Ingenito, E. P. (2003). On the Progressive Nature of Emphysema. *Am. J. Respir. Crit. Care Med.* 168, 516–521. doi:10.1164/rccm.200208-908PP

Sutton, M. A., Ke, X., Lessner, S. M., Goldbach, M., Yost, M., Zhao, F., et al. (2008). Strain Field Measurements on Mouse Carotid Arteries Using Microscopic Three-Dimensional Digital Image Correlation. *J. Biomed. Mater. Res.* 84A, 178–190. doi:10.1002/jbm.a.31268

Tawhai, M. H., and Bates, J. H. T. (2011). Multi-scale Lung Modeling. *J. Appl. Physiol.* 110, 1466–1472. doi:10.1152/japplphysiol.01289.2010

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. doi:10.1056/NEJM200005043421801

Tikhonov, A. N., Goncharsky, A., Stepanov, V. V., and Yagola, A. G. (2013). *Numerical Methods for the Solution of Ill-Posed Problems*. Springer Science & Business Media.

Toshima, M., Ohtani, Y., and Ohtani, O. (2004). Three-Dimensional Architecture of Elastin and Collagen Fiber Networks in the Human and Rat Lung. *Arch. Histology Cytol.* 67, 31–40. doi:10.1679/aohc.67.31

Vlahakis, N. E., Schroeder, M. A., Limper, A. H., and Hubmayr, R. D. (1999). Stretch Induces Cytokine Release by Alveolar Epithelial Cells *In Vitro*. *Am. J. Physiology-Lung Cell Mol. Physiol.* 277, L167–L173. doi:10.1152/ajplung.1999.277.1.L167

Wall, W. A., Wiechert, L., Comerford, A., and Rausch, S. (2010). Towards a Comprehensive Computational Model for the Respiratory System. *Int. J. Numer. Meth. Biomed. Engng.* 26, 807–827. doi:10.1002/cnm.1378

Werner, R., Ehrhardt, J., Schmidt, R., and Handels, H. (2009). Patient-Specific Finite Element Modeling of Respiratory Lung Motion Using 4D CT Image Data. *Med. Phys.* 36, 1500–1511. doi:10.1118/1.3101820

Zeng, Y. J., Yager, D., and Fung, Y. C. (1987). Measurement of the Mechanical Properties of the Human Lung Tissue. *J. Biomech. Eng.* 109, 169–174. doi:10.1115/1.3138661

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.

Edited by:

Jerome Noailly, Pompeu Fabra University, SpainReviewed by:

Abdelwahed Barkaoui, International University of Rabat, MoroccoBehzad Shariat, Université Claude Bernard Lyon 1, France

Thiranja Prasad Babarenda Gamage, The University of Auckland, New Zealand

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