All experiments conformed to the ethical principles and guidelines under protocols set forth and approved by the Boston University Institutional Animal Care and Use Committee (protocol number PROTO201900086). All animal procedures were compliant with ARRIVE guidelines. Mice were housed in ambient temperature and humidity and 12-h light–dark conditions under pathogen-free conditions at the Boston University Animal Science Center. No housing or handling exceptions were made for this study.

We used 11- to 23-week-old male and female mice for experimental procedures including healthy lung imaging and generating models of primary cancer, as previously described (Banerji et al., 2023). A breeding pair of transgenic B6.129(Cg)-Gt (ROSA)26Sortm4 (ACTB-tdTomato,-EGFP)Luo/J mice (JAX, 007676, Jackson Labs) (Muzumdar et al., 2007), referred to by the abbreviation “mTmG”, was initially purchased to breed a colony; that colony was the source of all animals for healthy lung and primary cancer experiments. For the present study, which examines two representative mice from this colony, the healthy mouse was 11 weeks old at the time of imaging. The urethane mouse, serving as our model of primary cancer, was 23 weeks at the start of urethane dosing and 54 weeks at the time of imaging. Between these ages, the murine lung’s volume does not change appreciably both in our experience and per development studies (Schulte et al., 2019).

We adapted a previously described protocol (Janker et al., 2018; Sozio et al., 2021) to induce primary lung cancer in mTmG mouse lungs using urethane (Sigma U2500). A stock solution of urethane was prepared at a working concentration of 200 mg/mL in PBS. Mice were dosed with the urethane solution at 1 mg/g body weight, twice weekly for 5 weeks by intraperitoneal (IP) injection. Mice were sacrificed and lungs harvested for imaging in the crystal ribcage after 6–12 months. The maximum tumor size permitted for the study was 1.5 mm in diameter. Mice were excluded from the study after presenting with labored breathing, hunched posture, or ruffled fur due to tumor progression.

The full development of the crystal ribcage platform is described in our previous work (Banerji et al., 2023). Briefly, microCT scans of C57BL/6 mouse chest cavity (courtesy the Hoffman group at the University of Iowa (Thiesse et al., 2010; Vasilescu et al., 2012; Kizhakke Puliyakote et al., 2016)) were segmented and refined to create the native ribcage geometry. In successive additive manufacturing and fabrication steps the ribcage model was converted into the crystal ribcage mold that was thermoformed over to create the polystyrene crystal ribcage. The internal surface was engineered to be hydrophilic to allow the lung to glide over its surface, as in the native ribcage. A six degree of freedom arm was included to rotate the crystal ribcage about any axis to image across the entire the distal lung surface using either a top-down or bottom-up configured microscope. Because lung volume changes significantly with age (Schulte et al., 2019), we have fabricated different, age-specific crystal ribcages to accommodate lungs of different sizes (Banerji et al., 2023).

Isolated mouse lungs were ventilated and perfused as previously described (Vanderpool and Chesler, 2011; Banerji et al., 2023). Briefly, the mouse trachea was cannulated and the lungs dynamically ventilated (Kent Physiosuite Mouse Ventilator, Kent Scientific). The lungs were perfused by cannulating the pulmonary artery and left atrium, and perfusing serum-free RPMI cell culture medium (Corning) through the lung vasculature. After cannulating the trachea and mouse heart, the lung–heart bloc, was excised and placed into the crystal ribcage for ex vivo microscopy under variable quasi-static positive air pressures.

As previously described (Banerji et al., 2023), ex vivo lungs, under quasi-static inflation conditions and within the crystal ribcage, were imaged using (i) an upright Nikon stereomicroscope with a 1x objective, and (ii) an upright Nikon CSU-X1 spinning-disk confocal microscope with 1x, 2x, 4x and 10x objectives, using NIS-Elements acquisition software and with the environmental temperature control set to 37°C.

Z-stacks of the diseased and healthy lungs were acquired on the confocal microscope using a 561 nm laser at 20–50 ms exposure (50–20 frames per second) per frame. Voxel sizes varied based on objective used, with XY resolution varying from 1–10 μm and Z step sizes varying from 2.5–12.5 μm. Total Z-stack acquisition time was on the order of 6–15 s for each positive-end expiratory pressure (PEEP) condition.

Before imaging, lungs were gradually recruited by slowly raising the intratracheal pressure to 18 cmH₂O, measured using custom sensors sensitive to 0.1 cmH₂O, using a water column. The pressure was then reduced in decrements of 1 cmH₂O down to 2 cmH₂O. To allow the lung to relax to its steady-state condition, each pressure was maintained for 1 min before imaging.

Due to significant variations in the lung’s intensity within a microCT volume, segmenting the organ by thresholding is unreliable. To segment the lung, we thus construct a naïve, Bayesian classifier—trained on a single, two-dimensional slice of the image along with its ground-truth class labels—to differentiate the lung class L from its surroundings M (Bishop, 2006). Let Ω ⊂ R 3 be the position vectors of pixels within the domain of the microCT image (Thiesse et al., 2010; Vasilescu et al., 2012; Kizhakke Puliyakote et al., 2016), so that x ⇀ ∈ Ω is the position of a given voxel. The image can then be expressed as the scalar field I x ⇀ : Ω → R . Let F x ⇀ : Ω → R n map from each of these position vectors to the n-dimensional feature vectors extracted from I x ⇀ : Ω → R . In the present study, each feature vector has ten components, φ i x ⇀ : Ω → R , each corresponding to a transformation of the image volume by a different neighborhood operation. These components are listed in Table 1.

Let Ω L ⊂ Ω be the subset of position vectors from a given slice of the image volume that have been manually labeled as belonging to the lung, and let Ω M ⊂ Ω be the remaining position vectors from the same slice. For the priors, we assume that the prior probability P L = Ω L Ω L + Ω M , which implies by the axiom of normalization that P M = 1 − P L . To estimate the likelihoods, P x ⇀   |   L and P x ⇀   |   M , we train Gaussian mixture models on the feature vectors F Ω L and F Ω M , respectively. Finally, from these definitions, we apply Bayes’ theorem to recover P L     x ⇀ , the posterior probability that a given pixel corresponds to lung tissue. P L     x ⇀ = P x ⇀     L P L P x ⇀     L P L + P x ⇀     M P M (1)

The lung segmentation S L x ⇀ : Ω → 0 , 1 is then defined as S L x ⇀ = 1 if P L     x ⇀ > 0.5 and S L x ⇀ = 0 otherwise. In the MATLAB implementation, these maps— I x ⇀ , F x ⇀ , S L x ⇀ , and P L     x ⇀ —are represented as matrices.

In contrast, because the intensity of bone tissue is much higher than that of other biological materials, the approach to segmenting the ribcage is simpler. Here, the ribcage segmentation S R x ⇀ : Ω → 0 , 1 is defined as S R x ⇀ = 1 if and only if the intensity I x ⇀ exceeds some constant, volume-dependent threshold. For both the lung segmentation and the ribcage segmentation, the resulting segmentation contains multiple connected components, corresponding to features like the scapula, humerus, cartilaginous tracheal rings, and tissue outside of the lung; extracting the largest connected components from these initial segmentations isolates the desired region of interest.

From the lung segmentation S L x ⇀ , we approximate the surface of the diaphragm as follows. First, we find z max x , y = ⁡ max     z     x , y , z ∈ Ω ⋀ S L x , y , z = 1 . We then filter the mapping z max x , y using a mode filter. Finally, we resample points from this surface using a thin-plate smoothing spline and save the point cloud to a text file (Hastie et al., 2009). This point cloud represents the geometry of the diaphragm.

Next, to construct a point-cloud approximation of the ribcage, we first find the geometric centroid of the lung, x ⇀ c , from S L x ⇀ using the following equation. x ⇀ c = 1 ∑ x ⇀ ∈ Ω S L x ⇀ ∑ x ⇀ ∈ Ω x ⇀   S L x ⇀ (2)

For each slice of the volume, we then project rays from the projection of the centroid, x ⇀ c , onto the given slice at a dense collection of angles from 0 to 2 π until each ray contacts a nonzero pixel on the interior of the ribcage segmentation S R x ⇀ . A thin-plate smoothing spline is then fit to these contact points to produce a surface approximating the ribcage, and a dense collection of points, P R , are sampled from this surface. Because the ribcage is open near the apex of the lung, these sampled points are artifactually peaked in the neighborhood of the apex. To correct this, the nodes near the apex are flattened by minimizing the following objective function. ∑ i k z i − Z i − g z 2 + α ∇ x y z i 2 + β ∇ x y 2 z i 2 + γ ∇ x y z i − Z i 2 + δ ∇ x y 2 z i − Z i 2 (3)

In this equation, z i represents the z-component of the ith node’s position vector after correction, Z i represents the same component before correction, k represents the stiffness of a virtual spring anchoring a point to its original height, and g z is a body force pulling these points toward the centroid of the ribcage. The remaining terms serve to regularize the optimization, penalizing the first and second derivatives of the height as well as changes in these derivatives. The sum is taken over the points near the apex of the lung. Finally, these point clouds of the diaphragm and ribcage are saved as text files for subsequent import into SolidWorks.

From these point clouds, we finally construct STEP (Standard for the Exchange of Product model data defined by ISO 10303 (Pratt, 2001)) representations of the diaphragm and the ribcage using the ScanTo3D feature in SolidWorks. By cutting the ribcage surface with the diaphragm surface, and filling the space enclosed between them, we recover a simplified model of the lung that is everywhere tangent to the ribcage. This approach guarantees a priori that, at the start of each simulation, the lung geometry and the ribcage geometry are in perfect contact, preventing errors in predicted strains and stresses that may arise due to mismatch between these geometries. The STEP representations of the ribcage and the lung are then exported from SolidWorks.

To perform finite-element simulations of the organ, these STEP geometries are now imported into Abaqus 2022 (Dassault Systèmes). The ribcage is taken to be a discrete, rigid part and is meshed with rigid, triangular elements. The lung is taken to be a deformable part and is meshed with C3D10 quadratic tetrahedral elements, which are chosen over C3D4 linear tetrahedral elements for their tendency to converge more quickly with coarser mesh resolutions. In simulations of the healthy lung, the lung mesh consists of 34,541 nodes and 21,945 elements; the ribcage mesh consists of 64,093 nodes and 127,697 elements. The simulation was performed on the Boston University Shared Computing Cluster hosted by the Massachusetts Green High-Performance Computing Center distributing the load over 12 processors with 4 GB of RAM per processor, each simulation completed within 7 h.

Based on prior studies and on the lung’s microstructure and constitutive behavior—which resembles a hyperelastic, tetrakaidekahedral foam whose walls are comprised of elastin, type-I collagen, and type-III collagen—the lung is modeled using the Ogden-Hill model (Berezvai and Kossa, 2017) of a hyperelastic foam (Vawter et al., 1979; Andrikakou et al., 2016). The general form of the strain-energy density function is thus taken to be U = ∑ i = 1 N 2 μ i / α i 2 λ 1 α i + λ 2 α i + λ 3 α i − 3 + 1 β i J − α i β i − 1 , (4) where α i is a dimensionless material parameter determining the nonlinear behavior of the stress-strain relation, β i is a dimensionless material parameter given by the Poisson’s ratio as β i = ν i / 1 − 2 ν i , μ i is a material parameter with units of stress determining the shear modulus during small strains from the reference configuration, λ i is the i t h principal stretch, J is the determinant of the deformation gradient, and N is the number of terms in the model. For simplicity, we assume that the strain-energy density function consists of only one term, reducing the general equation to U = 2 μ / α 2 λ 1 α + λ 2 α + λ 3 α − 3 + 1 β J − α β − 1 . (5)

This strain-energy density function, and consequently the pressure, is linear in the parameter μ  and nonlinear in the parameter  α . Therefore, if we simulate the distensions for parameters ( μ , α , then the pressures required to produce the same distensions for any other μ * , α are a simple rescaling of those for ( μ , α . In practice, therefore, it is unnecessary to simulate ventilation for the same α but different μ to determine the pressure-distension curve.

The boundary conditions (Supplementary Figure S2) include immobilization of the ribcage, frictionless sliding contact between the lung’s upper surface and the ribcage, and negative pressure on the boundary of the lung. Although experiments involve positive-pressure ventilation, the simulation involves applying negative pressure to the external surface of the lung; the explanation for this apparent discrepancy is that the governing equations are symmetric under mutual inversion of the pressure’s sign and the surface normal’s direction, implying that the model is equally applicable to either mode of ventilation (Shi et al., 2022). Since the parts have been designed a priori to be tangent everywhere, initial contact between the surfaces is easy to establish.

While previous studies (Tawhai et al., 2009; Shi et al., 2022) have shown that gravity significantly influences the mechanics of the human lung, our model neglects the influence of gravity due to its smaller role in the mouse. Consider the conservation of linear momentum under conditions of static equilibrium, which has been rendered dimensionless (Munson et al., 2012) by factoring out the lung density ρ lung , gravitational acceleration g , lung height H lung , and transpulmonary pressure P tp . P tp ρ lung g H lung ∇ * ⋅ P * + B * = 0 (6)

Here, ∇ * , P * , and B * represent the dimensionless divergence operator, the dimensionless tissue stress, and the dimensionless gravitational body forces, respectively. The tissue stress is known to be approximately equal to the transpulmonary pressure (Mead et al., 1970). Consequently, the dimensionless number P tp ρ lung g H lung characterizes the magnitude of the tissue stress divergence relative to the magnitude of gravity. In the human lung, this dimensionless quantity remains below 1 for transpulmonary pressures up to 20 cmH₂O, indicating the significant role of gravity in governing its mechanics. In the mouse lung, however, this dimensionless quantity is equal to 1 when the transpulmonary pressure is 1 cmH₂O, but decreases linearly as the transpulmonary pressure increases; once the transpulmonary pressure reaches 10 cmH₂O, this dimensionless number increases to 10, indicating that tensile forces within the tissue greatly exceed gravitational body forces. Based on this reasoning, we posit that it is reasonable to neglect gravity when modeling the murine lung, with the approximation improving at higher pressures. Additional reasoning is discussed in the results.

Finally, it is important to note that the lung exhibits hysteresis, with its inflation characterized by one strain-energy density function and its deflation characterized by another (Vawter et al., 1979). In this study, we elect to model the lung’s behavior during quasistatic deflation, so that our measurements used to calibrate the model are taken from states of higher pressure to states of lower pressure in near-equilibrium. The same approach can easily be repeated to recover a model of the lung’s behavior during inflation.

To simulate the cancerous lung, the same process is repeated, but rather than being constant, the material field μ X ⇀ is defined as μ X ⇀ = A e − X ⇀ − X ⇀ tumor n / σ n + B , (7) where X ⇀ tumor is the centroid of the spherical tumor in material coordinates, σ determines its width, n controls the shape of the decay in magnitude with distance, and the coefficients A and B determine the stiffness in the near and far fields. To illustrate, we note that μ X → approaches A + B as X → approaches X → tumor , and that μ X → approaches B as X → − X → tumor approaches infinity. Consequently, A + B is the stiffness at the tumor’s centroid, while B is the stiffness far from the tumor. This material equation possesses three attractive properties: (i) it is spherically symmetric, corresponding to our interpretation that the equation represents a spherical tumor, (ii) stiffness decays with distance from the centroid, reflecting our intent that the tumor is stiffer than its surroundings, and (iii) it is differentiable and smooth, making it easier to work with during numerical computations. When n = 2 , the field is a multidimensional normal distribution, and as n → ∞ , the field approaches an indicator function for a ball. In our studies, we chose n = 4 .

The ribcage mesh is the same as before, while the lung mesh now consists of 12,679 nodes and 7,766 elements. Adaptive mesh refinement, which was necessary for convergence in the neighborhood of the tumor, yielded a mesh with fewer elements relative to the simulation of the healthy lung. The simulation was performed on the Boston University Shared Computing Cluster hosted by the Massachusetts Green High-Performance Computing Center; distributing the load over 12 processors with 4 GB of RAM per processor, the simulation completed within 4 h.

The finite-element simulation is repeated for different material parameters within a neighborhood of the values yielding optimal agreement between empirical and in silico outcomes. When the strain-energy density function is restricted to a single term, the parameter μ represents the material’s shear modulus during small strains, as can be shown by a Taylor expansion of the strain-energy density function about λ = 1 . From the stress-strain curve in the linear regime, when the strain increases linearly with distending pressure, we can estimate the Young’s modulus as the ratio of the increment in strain to the increment in pressure. From Figure 2C, we estimate that incrementing the pressure from 0 kPa to 0.5 kPa yields an increment in strain of 0.2. Following this reasoning, we find that E ≈ 0.5  kPa / 0.2 = 2.5  kPa . Assuming a Poisson’s ratio of ν = 0.2 (Concha et al., 2018), this value for the Young’s modulus corresponds to a shear modulus of μ = E / 2 1 + ν ≈ 1   k P a .

Based on these initial estimates, we perform simulations for μ = 1  kPa and the dimensionless parameter α ranging from six to 16. The outcomes of these simulations are then used to construct a calibration surface P λ H , μ , α , where λ H is the stretch of the finite-element model along the vertical axis, μ and α are the material parameters mentioned earlier, and P is the transpulmonary pressure required to distend the finite-element model to this degree of stretch. From the observed distensions of the lung in the crystal ribcage (Figure 1B) for pressures ranging from 2 to 18 cmH₂O, λ o b s , P o b s , we then determine the maximum likelihood assignments for μ and α by minimizing the following objective function (Bishop, 2006). ∑ λ o b s , P o b s P λ o b s , μ , α − P o b s 2 (8)

This optimization procedure is not computationally intensive; running on a personal laptop with 16 GB of RAM and a typical CPU, for example, the process completes within seconds. The outcome is a single, point estimate of the coefficients that characterize the organ-scale behavior. It should be noted that, because the model is linear in the parameter μ , we can in practice simplify the problem of determining the optimal parameters by transforming the objective function as follows. ∑ λ o b s , P o b s μ P λ o b s , 1 , α − P o b s 2 (9)

This implies that, rather than sampling the function P λ o b s , μ , α , we only need to sample the function P λ o b s , 1 , α in our finite-element simulations. In contrast, because the model is nonlinear in α , an analogous transformation is not possible for α ; this is why we limit the constitutive model to a single term.

To measure the displacements caused by a change in distending pressure applied to the lung at cellular resolution, we leverage deformable image registration (Heinrich et al., 2012a) (Figure 1C). To prepare the images for registration, we perform an optimization to autonomously correct the bulk rotation of the material due to the natural curvature of the crystal ribcage away from the imaging plane. Because the image-registration algorithm does not, in practice, yield good estimates for the displacements along the shallow depth axis, we project the rotationally corrected images along this axis.

After preprocessing the images, we invoke the image-registration algorithm. This algorithm, which has been adapted from earlier literature (Heinrich et al., 2012a), formulates the inverse-elasticity problem as an inference problem on a Markov Random Field (MRF) to automatically determine the displacements necessary to match images of the lung at two different pressures. The observable variables of the Markov Random Field are modality-independent neighborhood descriptors (MIND) extracted from the image (Heinrich et al., 2012b), while the latent variables are the displacements necessary to minimize the sum of squared differences in these descriptors across the two images; edges between the latent variables of neighboring observables represent the constraint that displacements vary smoothly throughout the domain. Let O and T be the original and deformed images, respectively. Furthermore, let W O , u be the image produced by warping the image O using the displacements u . Lastly, let the function D map an image to its MIND representation. Then the registration algorithm finds the displacements that minimize the objective function ∑ x ⇀ ∈ Ω D W O , u −   D T 2 , (10) where the sum is taken over the points in the image domain, Ω . From these displacements and the material properties of the calibrated finite-element model, we next find the distribution of stiffnesses throughout the domain, as described in the next section.

To determine the relative stiffnesses of the material throughout the image domain, we leverage a Python implementation of the Adjoint-Weighted Equation (AWE) formulation of the inverse-elasticity problem (Barbone et al., 2010). Let X ⇀ ∈ Ω 0 ⊂ R 2 be the coordinates of material points of a deformable body in the reference configuration, and let x ⇀ ∈ Ω 1 ⊂ R 2 be the coordinates of the same material points after some deformation. Because soft tissues deform continuously, there exists a continuous function ψ : R 2 → R 2 such that x ⇀ = ψ X ⇀ . Given the registered displacements, u X ⇀ = ψ X ⇀ − X ⇀ , along with the definition of the deformation gradient, F X ⇀ = ∇ X ⇀ ψ X ⇀ , we can recover the deformation gradient tensor in terms of X ⇀ (Holzapfel, 2000). F X ⇀ = I + ∇ X ⇀ u X ⇀ (11)

From the deformation gradient, we next determine the Cauchy-Green strain tensor field, C X ⇀ , as follows. C X ⇀ = F X ⇀ T F X ⇀ (12)

By the spectral theorem and the manifest symmetry of C X ⇀ , we subsequently determine the eigenvalues, λ 1 X ⇀ 2 and λ 2 X ⇀ 2 , of this tensor.

The areal strain, which is a scalar field representing the fractional change in the material’s area relative to its value in the reference configuration, is then given by ε A X ⇀ = λ 1 X ⇀ λ 2 X ⇀ − 1 . (13)

The whole organ’s stress-strain behavior is well-described by a hyperelastic material law, where the nominal stress is the derivative of the strain-energy function with respect to the principal nominal stretches. If we assume that this model holds at all length scales, down to the cellular scale and for both healthy and diseased tissue, then we can use this same constitutive model when solving the inverse problem. Therefore, from the same strain-energy function introduced earlier, we recover the principal components of the first Piola-Kirchhoff stress tensor as follows. P i = ∂ W ∂ λ i = μ X ⇀ 2 α λ i λ i α − J − α β (14) where β = ν 1 − 2 ν . In evaluating the above expression, we need to compute the Jacobian determinant J = λ 1 λ 2 λ 3 , but image registration only yields the stretches tangential to the lung’s surface. Consequently, we approximate J as J ≈ λ 1 λ 2 3 / 2 . Finite-element simulations of the whole organ indicate that this approximation generally holds within 10%–15% error (Supplementary Figure S2).

For a hyperelastic material, it is well-known that the eigenvectors of the stress are aligned with the eigenvectors of the strain. Therefore, if v ⇀ i are the principal directions of the right Cauchy-Green strain tensor, the first Piola-Kirchhoff stress tensor becomes P = P 1 v ⇀ 1 ⊗ v ⇀ 1 + P 2 v ⇀ 2 ⊗ v ⇀ 2 . (15)

We can simplify the above expression by defining the tensor A X ⇀ = P X ⇀ / μ X ⇀ , leading to the following simplified form. P X ⇀ = μ X ⇀ A X ⇀ (16)

Here, μ X ⇀ is unknown while A X ⇀ is completely determined by the displacements. At static equilibrium and in the absence of body forces, the conservation of linear momentum requires that the divergence of the first Piola-Kirchhoff stress, P X ⇀ , with respect to the material coordinates, X ⇀ , vanishes. ∇ X ⇀ ⋅ P X ⇀ = 0 (17)

From our earlier simplified form for P X → , the condition of static equilibrium becomes ∇ X ⇀ ⋅ μ X ⇀ P X ⇀ = 0 . (18)

In the AWE formulation of the inverse elasticity problem for a linearly elastic material, we seek a variational solution for μ X → to the above differential equation. This yields a map μ X ⇀ of relative stiffnesses, whose scale is determined by the specified-mean boundary condition (Albocher et al., 2009). After solving the inverse problem for μ X ⇀ , we rescale μ X ⇀ so that the mean stiffness throughout the domain matches the organ-scale stiffness determined by calibrating the finite-element model.

From Holzapfel’s equation (6.180) (Holzapfel, 2000), which is reproduced below as Eq. 19, along with the strain-energy function mentioned earlier, we determine the components of the elasticity tensor, C , at each state of deformation for each position within the domain. In the following equation, S a are the principal values of the second Piola-Kirchoff stress tensor, λ a are the principal stretches, and N ^ a are the principal directions of the deformation. C = ∑ a , b = 1 3 1 λ b ∂ S a ∂ λ b N ^ a ⊗ N ^ a ⊗ N ^ b ⊗ N ^ b + ∑ a , b = 1 , a ≠ b 3 S b − S a λ b 2 − λ a 2 N ^ a ⊗ N ^ b ⊗ N ^ a ⊗ N ^ b + N ^ a ⊗ N ^ b ⊗ N ^ b ⊗ N ^ a (19)

Finally, we perform an iterative optimization to determine the Young’s modulus and Poisson’s ratio fields that optimally approximate the components this tensor. This Young’s modulus is what we ultimately report as the lung’s stiffness (Figure 1D).