1 Introduction

Around 334 million individuals worldwide suffer from asthma and it is estimated that over 250,000 of these people die prematurely each year as a result (Forum of International Respiratory Societies 2017). Asthma is the most prominent chronic disease amongst youths, affecting over 14% of children globally (Pearce et al. 2007), yet despite its rising prevalence, the cause and onset of asthma remains unknown. Understanding asthma is of vital importance.

Asthma is characterised by inflammation, airway hyperresponsiveness and remodelling (Brightling et al. 2012; Berair et al. 2013). Airway hyperresponsiveness refers to excessive bronchoconstriction (narrowing of the airway) due to rapid contraction of airway smooth muscle (ASM) in response to a relatively low dose of contractile agonist (King et al. 1999; West et al. 2013). Chronic inflammation causes swelling of the airway tissue, narrowing the airway (León 2017), and resulting in overall restricted pulmonary function. The persistent structural changes due to inflammatory injury repair, airway thickening and scarring constitute airway remodelling (Bossé et al. 2008; Al Alawi et al. 2014). Until recently, airway remodelling has been predominantly attributed to chronic inflammation (Saglani and Lloyd 2015). Current experimental evidence, however, suggests that bronchoconstriction induced airway narrowing may play a key role in promoting remodelling (Grainge et al. 2011) via activation of the regulatory cytokine, transforming growth factor \(\beta \) (TGF-\(\beta \)) (Wipff et al. 2007; Buscemi et al. 2011; Tatler et al. 2011). In addition, TGF-\(\beta \) has been shown to act as a contractile agonist (Ojiaku et al. 2018). Precision-cut lung-slice (PCLS) stretching experiments are a widely-used ex vivo assay (see, e.g., Sanderson (2011) and Tan and Sanderson (2014)) for studying agonist-driven bronchoconstriction and more recently how this links to mechanical activation of TGF-\(\beta \). In particular, Tatler (2016) showed that stretching of PCLS increases TGF-\(\beta \) activation. However, results from such studies are difficult to interpret without knowledge of the underlying tissue deformation and stress state.

Airway mechanical behaviour is dominated by ASM and collagen-rich extracellular matrix (ECM). ASM is arranged in fibrous bundles that are oriented helically within the airway wall; this arrangement is thought to enhance the bronchocontractile ability of the smaller airways (Amrani and Panettieri 2003; Ijpma et al. 2017). The ECM is a delicate mesh of deposited connective tissue and fibrous proteins (Hinz 2015; Cheng et al. 2016), surrounding ASM cell bundles. In the undeformed state, collagen exhibits a ‘crimped’ structure (Kadler et al. 1996); under strain, these structures straighten to bear load. Varition in their natural lengths mean that ECM is ‘recruited’ successively as strain is increased, imparting a strain-stiffening behaviour to the airway (Wells 2013).

Early mathematical models of the intact airways, accounting for tension generated by ASM contraction and mechanical properties of the airway wall (e.g. Latourelle et al. (2002), Affonce and Lutchen (2006), and Ma and Lutchen (2006)), were based on empirical stress-strain relationships for the whole airway (Lambert et al. 1982, 1993, 1994; Lambert and Wilson 2005) and the Laplace thin-airway wall approximation (Anafi and Wilson 2001). With these models, it is not possible to determine tissue stresses within the airway wall nor to separate ASM and ECM contributions to the mechanics. Similarly, early models of airways embedded in parenchyma mimicking PCLS experiments assumed the thin-wall Laplace approximation (Bates and Lauzon 2007; Khan et al. 2010). Brook et al. (2010) extended these early models to account for multiple airway constituents assuming a finite airway wall thickness for both the intact airway and PCLS models under a plane strain and plane stress approximation respectively. While these allowed for tissue stresses to be determined within the airway wall, the linear elastic framework used meant that predictions were only qualitatively useful. Breen et al. (2012) considered a finite element finite-elasticity model of the PCLS but neglected airway wall thickness and was only concerned with stresses in the lung parenchyma. Following arterial and cardiovascular mechanics (Gasser et al. 2006; Ateshian 2007; Holzapfel and Ogden 2010; Hill et al. 2012), a nonlinear elastic single-phase fibre-reinforced airway model assuming finite airway wall thickness under a plane strain approximation was developed by Hiorns et al. (2014). Examples of models that account for mucosal growth, buckling and folding are: Wiggs et al. (1997), Moulton and Goriely (2011), Li et al. (2011), Eskandari et al. (2015). More generally, the mechanics of growth in thin biological membranes is described by Kroon and Holzapfel (2008) and Rausch and Kuhl (2014). Approximate solutions for axisymmetric stretching of thin elastic membranes with traction-free surfaces have previously been determined for isotropic materials (Wong and Shield 1969; Yang 1967) but do not account for anisotropy or active contraction. Others consider finite deformations of incompressible rubber membranes with a central solid inclusion and under uniform pressure (Jianbing et al. 2015). Finally, models of agonist driven feedback that are focused on growth and remodelling are presented by Chernyavsky et al. (2014), Aparício et al. (2016), and Hill et al. (2018). However, none of these previous models are suitable descriptions for finite deformation of a thin slice in which the tissue stresses may be determined.

In this study we develop a nonlinear fibre-reinforced biomechanical model of PCLS stretching experiments accounting for ASM contraction in response to agonist exposure and ECM strain-stiffening. Through numerical simulation, we quantify the mechanical stress experienced by the airway wall constituents in response to cyclic stretching that consequently activates TGF-\(\beta \). Additionally, we assess the applicability of two simplifying limits of the model; namely, a one-dimensional membrane representation of the PCLS and an asymptotic reduction in the thin-slice limit (that nevertheless retains a description of axial deformation).

2 A biomechanical model of PCLS

The PCLS is a well-established experimental preparation for studying airway reactivity, and corresponding biomechanical response (see, e.g., Wang et al. (2008) and Tan and Sanderson (2014)). The key advantage of the PCLS is that vital functional interactions between airways, arterioles, and veins are preserved within the alveolar parenchyma (Sanderson 2011). PCLS are obtained by inflating human lung tissue with liquid agarose, which is allowed to set and solidify before finely slicing. Stretching of the PCLS is effected by adhering it to a deformable membrane, to which a stretch is applied (Fig. 1). In the experiments of Tatler (2016), that form our primary motivation, stretch is applied cyclically, in the form of a sine wave with a 15% amplitude and 0.3Hz frequency, for 24 hours (Tatler 2016). Stretch is applied with a 5%, 10% and 20% amplitude thereafter.

We represent a single airway within the PCLS as a cylinder, whose constituents are modelled via constrained mixture theory (Truesdell and Toupin 1960; Bowen 1976; Truesdell and Noll 1992; Ateshian 2017). The formulation for obtaining the constitutive mechanical relation for this type of material is given in detail by Humphrey and Rajagopal (2002), Ateshian (2007) and Ateshian and Ricken (2010). Specifically, we consider a saturated multiphase mixture of an active contractile ASM component and a passive ECM component, each modelled as a nonlinear, incompressible, fibre-reinforced hyperelastic material (Holzapfel 2000), with associated volume fractions

$$\begin{aligned} \Phi _c&= \frac{c^*}{{\bar{c}^*}},&\Phi _m&=\frac{m^*}{{\bar{m}^*}}, \end{aligned}$$
(1)

respectively, where \(c^*\) and \(m^*\) are the apparent densities and \(\bar{c}^*\) and \(\bar{m}^*\) denote the true densities of ASM and ECM, respectively. The assumption of intrinsic incompressibility and tissue saturation demands

$$\begin{aligned} \sum _{i} \Phi _i&= 1,&i&\in \{c,m\}. \end{aligned}$$
(2)

ECM strain-stiffening occurs in the direction of the collagen fibre orientation and accounts for the recruitment of collagen fibres (from a crimped to uncrimped configuration) when stretched (Hiorns et al. 2014). Contractile force generation is assumed to occur in the direction of the ASM bundle orientation and occurs in response to an exogenous agonist and/or active TGF-\(\beta \) signalling pathways. Since the duration of the PCLS experiment of interest is significantly less than that of ASM growth or proliferation and ECM deposition, we assume that \(\Phi _c\) and \(\Phi _m\) are constant. In addition, for simplicity we neglect the time-dependent feedback between tissue strain and TGF-\(\beta \) activation. Furthermore, for simplicity, tissue porosity and constituent volume fraction changes that arise through large deformations are not considered in this study.

Fig. 1
figure 1

(i) Axisymmetric cyclic stretching of PCLS via the BioFlex method (Tatler 2016). A PCLS is adhered to a circular deformable rubber membrane and then an axisymmetric cyclic stretch is applied to the membrane (via a vacuum) in order to stretch the attached PCLS. (ii) Representative image of a murine PCLS adhered to the flexible membrane of a BioFlex culture well plate prior to stretching regime. (iii) Representative image of an intact airway within lung tissue of a murine PCLS post adherence to a BioFlex culture well (x20 magnification) showing airway and airspaces. Lung tissue was stained for alpha-smooth muscle actin (brown stain) to highlight smooth muscle cells. (iv) Dimensional undeformed reference configuration (left) and deformed configuration (right) to illustrate the geometry of an airway modelled within the PCLS. Dotted lines indicate the circumferential fibres representing the active ASM and passive ECM components

Full size image

2.1 Geometry and constitutive assumptions

Following the traditional continuum mechanics approach (Truesdell and Toupin 1960; Truesdell and Noll 1992; Holzapfel 2000) we assume that a common unstressed and unstrained reference configuration applies to each constituent in the airway, in which Lagrangian cylindrical coordinates \(\left( R^*,\Theta ,Z^* \right) \) describe the airway geometry:

$$\begin{aligned} R^*_{\text {in}} \le R^* \le R^*_{\text {out}}, \qquad 0 \le \Theta \le 2\pi , \qquad -\tfrac{h^*}{2} \le Z^* \le \tfrac{h^*}{2}, \end{aligned}$$
(3)

and wherein asterisks denote dimensional quantities, \(R^*_{\text {in}}\) and \(R^*_{\text {out}}\) denotes the inner and outer radius and \(\pm \frac{h^*}{2}\) denotes the upper and lower surfaces of the undeformed airway, respectively (see Fig. 1 (iv)).

Imposed stretching (herein we use stretch to imply elongation) or contraction of the ASM causes a deformation described by the deformed configuration \((r^*, \theta , z^*)\). For simplicity, we consider an axisymmetric radial airway stretch, and further assume there is no torsion, so that the deformation is described by

$$\begin{aligned} r^* = {r}^*(R^*,Z^*), \qquad \theta = \Theta , \qquad z^* = {z}^*(R^*,Z^*). \end{aligned}$$
(4)

The constituents within the tissue are constrained and therefore also deform axisymetrically according to (4).

The deformation gradient tensor for the tissue and each constrained constituent within is defined by \(\varvec{\mathbf{{F}}} \equiv \pmb {\nabla }{\mathbf {x}}\) in the reference configuration and in cylindrical polars is given by

$$\begin{aligned} \varvec{\mathbf{{F}}}= \begin{pmatrix} \frac{\partial r^*}{\partial R^*} &{}\quad 0 &{}\quad \frac{\partial r^*}{\partial Z^*} \\ 0 &{}\quad \frac{r^*}{R^*} &{}\quad 0 \\ \frac{\partial z^*}{\partial R^*} &{}\quad 0 &{}\quad \frac{\partial z^*}{\partial Z^*} \end{pmatrix}. \end{aligned}$$
(5)

The left and right Cauchy–Green deformation tensors are defined by \(\mathbf{B}=\mathbf{F}{} \mathbf{F}^{T}\) and \(\mathbf{C}=\mathbf{F}^{T}{} \mathbf{F}\), respectively. Incompressibility of the tissue is enforced by demanding (Ogden 2003)

$$\begin{aligned} \mathrm {det}[\mathbf{F}] = 1. \end{aligned}$$
(6)

Mechanical anisotropy is imparted to the airway via strain-stiffening of collagen fibres and contractile force generation of ASM bundles (Ogden 2003), as summarised above. The ASM and ECM constituents within the tissue are associated with a set of helically orientated fibers (see Fig. 1, Ijpma et al. (2017)). For simplicity, however, we describe this as a single set of fibres orientated circumferentially with undeformed direction denoted \({\mathbf {G}}\). In the deformed configuration the fibres have direction \({\mathbf {g}}=\mathbf{F}{\mathbf {G}}\).

Under the above assumptions, the constitutive mechanical law for the airway wall is obtained following the additive approach of Ambrosi and Pezzuto (2012) by introducing an active component, \(\Psi ^*_{\text {act}}\), to the passive isotropic, \(\Psi ^*_{\text {iso}}\), and anisotropic, \(\Psi ^*_{\text {ani}}\), components of the strain-energy function, \(\Psi ^*\). For each constituent within the tissue, we follow Holzapfel et al. (2000) and define the strain-energy for the ASM and ECM within the tissue as

$$\begin{aligned} {\Psi _i}^*(I_1,I_4) = {\Psi _i}^*_{\text {iso}}(I_1) + {\Psi _i}^*_{\text {ani}}(I_4) + {\Psi _i}^*_{\text {act}}(I_4), \qquad i \in \{c, m\}, \end{aligned}$$
(7)

wherein (here and throughout) the subscripts \(i\in \{c,m\}\) denote variables associated with each phase. In (7), \(I_1\) and \(I_4\) denote the first and fourth principle invariants of the left Cauchy-Green deformation tensor, \(\mathbf{B}\), and are defined by

$$\begin{aligned} I_1 = \mathrm {tr}[\pmb {{\mathsf {B}}}], \qquad I_4 = {\mathbf {g}}\cdot {\mathbf {g}}. \end{aligned}$$
(8)

We assume that the isotropic response of the tissue is described via a Neo-Hook-ean constitutive law, with passive isotropic stiffness \({\mu ^*_c}\) and \({\mu ^*_m}\) for the ASM and ECM, respectively. It is assumed that collagen fibres within the ECM do not store strain-energy when the airway is not inflated, i.e. at low transmural pressures, hence following Holzapfel et al. (2000) we associate an isotropic part of the ECM strain-energy function to the mechanical response of the non-collagenous matrix material. At high transmural pressures the resistance of the tissue to stretch is almost entirely due to anisotropic collagen fibre recruitment within the ECM. To account for strain-stiffening, as in Hiorns et al. (2014), we employ the anisotropic model of Holzapfel et al. (2000) with the addition of a Heaviside function so that the collagen fibres are only recruited when stretched. This approach provides a convenient way of abolishing the contribution of fibres in compression (i.e. when \(I_4 \le 1\)) and, as discussed in Holzapfel et al. (2004), satisfies the relevant necessary convexity constraints on such a strain-energy function. There is no active force contribution from the ECM; however, we include an active component in the ASM strain-energy function. The form of the active component in the Cauchy stress tensor (denoted \(\pmb {\sigma }^*\) and defined below) follows the general form used in Ambrosi and Pezzuto (2012),

$$\begin{aligned} {\pmb {\sigma }^*_c}_{\text {act}} = {\alpha }^*({\mathbf {g}}\otimes {\mathbf {g}}), \end{aligned}$$
(9)

where \({\mathbf {g}}\) denotes the direction of the deformed fibres and \(\alpha ^*\) is the active contractile force density (force per unit area) generated by the ASM. In reality the contractile force density will vary with both time-dependent stretch and/or time-dependent changes in exogenously applied contractile agonist. In this study, however, our focus is purely on examining the effect of a contractile force on the approximations made in steady state, and therefore, for simplicity we assume that \(\alpha ^*\) is constant.

In view of the above, the strain-energy functions for the ASM and ECM components are given by

$$\begin{aligned} {\Psi ^*_c}(I_1,I_4)&= \frac{{\mu ^*_c}}{2}(I_1 - 3) + \frac{{\alpha }^*}{2}I_4, \end{aligned}$$
(10a)
$$\begin{aligned} {\Psi ^*_m}(I_1,I_4)&= \frac{\mu ^*_m}{2}(I_1 - 3) + \frac{\omega ^*}{2\zeta } \mathrm {H}(I_4 -1)\left( \exp \left( \zeta (I_4 -1)^2\right) - 1 \right) . \end{aligned}$$
(10b)

Here, \(\omega ^* > 0\) is a constant parameter defining the passive anisotropic stiffness and accounts for the density of the fibres in the matrix and \(\zeta > 0\) is a dimensionless constant parameter defining the nonlinear increase in stiffness of the fibres as they deform (Hiorns et al. 2014). Differentiating (10) with respect to the invariants \(I_1\) and \(I_4\), respectively, we have

$$\begin{aligned} {\psi _c}^*_1&= \frac{{\mu _c}^*}{2}, \end{aligned}$$
(11a)
$$\begin{aligned} {\psi _c}^*_4&= \frac{{\alpha }^*}{2}, \end{aligned}$$
(11b)
$$\begin{aligned} {\psi _m}^*_1&= \frac{{\mu _m}^*}{2}, \end{aligned}$$
(11c)
$$\begin{aligned} {\psi _m}^*_4&= \omega ^*(I_4 - 1)\mathrm {H}(I_4 -1) \exp \left( \zeta (I_4-1)^2\right) , \end{aligned}$$
(11d)

where \({\psi _i}^*_j = \partial \Psi ^*_i / \partial I_j\) with \(i \in \{c,m\}\) and \(j \in \{1,4\}\).

The Cauchy stress tensor for each constituent is defined by,

$$\begin{aligned} \pmb {\sigma }^*_{i} = -{\mathscr {P}}^* \mathbf{I} + 2{\psi _i}^*_1\mathbf{B} + 2{\psi _i}^*_4 {\mathbf {g}}\otimes {\mathbf {g}}, \qquad i \in \{c,m\} . \end{aligned}$$
(12)

Here, \(\mathbf{I}\) is the identity matrix, the pressure \({\mathscr {P}}^*\) is a Lagrange multiplier included to enforce tissue incompressibility, and \({\psi _i}^{*}_{j}\) with \(i \in \{c,m\}\) and \(j \in \{1,4\}\) denote the derivatives of the strain-energy functions defined in (11).

Using (10), the strain-energy function for the whole tissue, \(\Psi ^*\), is

$$\begin{aligned} {\Psi }^*(I_1,I_4) = \Phi _c{\Psi _c}^*(I_1,I_4) + \Phi _m{\Psi _m}^*(I_1,I_4). \end{aligned}$$
(13)

and similarly, using (11), the strain-energy function derivatives for the whole tissue are,

$$\begin{aligned} \psi _1^*&= \Phi _c{\psi _c^*}_1 + \Phi _m{\psi _m^*}_1, \end{aligned}$$
(14a)
$$\begin{aligned} \psi _4^*&= \Phi _c{\psi _c^*}_4 + \Phi _m{\psi _m^*}_4. \end{aligned}$$
(14b)

The Cauchy stress tensor for the whole tissue, \(\pmb {\sigma }^*,\) is defined

$$\begin{aligned} \pmb {\sigma }^* = -{\mathscr {P}}^*\mathbf{I} + 2\psi _1^*\mathbf{B} + 2\psi _4^* {\mathbf {g}}\otimes {\mathbf {g}}, \end{aligned}$$
(15)

Equivalently, (15) may be obtained via the weighted sum of the Cauchy stress components for each constituent, (12), such that \(\pmb {\sigma }^* = \sum _i \Phi _i{\pmb {\sigma }_i}^*\).

2.2 Governing equations and boundary conditions

In mechanical equilibrium, and assuming there are no body forces on the tissue, the balance of linear momentum requires

$$\begin{aligned} \nabla \cdot \pmb {\sigma }^* = {\mathbf {0}}, \end{aligned}$$
(16)

subject to the following boundary conditions.

At the outer radius, we enforce a displacement boundary condition,

$$\begin{aligned} r^*\left( R^*_{\text {out}},Z^*\right) = r_{\text {dis}}^*, \end{aligned}$$
(17)

to mimic the axisymmetric stretch imposed on the PCLS via the BioFlex method (Fig. 1). For simplicity, time-dependent loading is not considered in this study and we assume that \(r_{\text {dis}}^*\) is a constant.

The upper, lower and inner surfaces of the tissue are traction-free such that

$$\begin{aligned} \pmb {\sigma }^*\left( R^*,\tfrac{h^*}{2}\right) {\mathbf {n}}^*_{\text {up}}&= {\mathbf {0}},\end{aligned}$$
(18a)
$$\begin{aligned} \pmb {\sigma }^*\left( R^*,-\tfrac{h^*}{2}\right) {\mathbf {n}}^*_{\text {low}}&= {\mathbf {0}},\end{aligned}$$
(18b)
$$\begin{aligned} \pmb {\sigma }^*\left( R_{\text {in}}^*,Z^*\right) {\mathbf {n}}^*_{\text {in}}&= {\mathbf {0}}, \end{aligned}$$
(18c)

where the unit normals to the upper, \({\mathbf {n}}^*_{\text {up}}\), lower, \({\mathbf {n}}^*_{\text {low}}\), and inner, \({\mathbf {n}}^*_{\text {in}}\), surfaces in the reference configuration are given by:

$$\begin{aligned} {\mathbf {n}}^*_{\text {up}, \text {low}}&= \left( \left( -\frac{\partial z^*}{\partial R^*}\right) ^2 + \left( \frac{\partial r^*}{\partial R^*}\right) ^2\right) ^{-\tfrac{1}{2}} \cdot \left( -\frac{\partial z^*}{\partial R^*} , 0,\frac{\partial r^*}{\partial R^*}\right) , \qquad \text {at~} Z^*=\pm \tfrac{h^*}{2};\end{aligned}$$
(19a)
$$\begin{aligned} {\mathbf {n}}^*_{\text {in}}&= \left( \left( -\frac{\partial z^*}{\partial Z^*}\right) ^2 + \left( \frac{\partial r^*}{\partial Z^*} \right) ^2\right) ^{-\tfrac{1}{2}} \cdot \left( -\frac{\partial z^*}{\partial Z^*}, 0,\frac{\partial r^*}{\partial Z^*} \right) , \qquad \text {at~} R^* = R^*_{\text {in}}. \end{aligned}$$
(19b)

2.2.1 Non-dimensionalisation

We non-dimensionalise the governing equations by introducing the following scalings

$$\begin{aligned} (r,R) = \frac{(r^*,R^*)}{R^*_{\text {out}}}, \qquad (z,Z) = \frac{(z^*,Z^*)}{h^*}, \qquad ({\mathscr {P}}, \Psi _i) = \frac{({\mathscr {P}}^*, \Psi _i^*)}{\mu ^*_c}, \end{aligned}$$
(20)

so that the dimensionless undeformed reference configuration is given by

$$\begin{aligned} R_{\text {in}} \le R \le 1, \qquad 0 \le \Theta \le 2\pi , \qquad -\tfrac{1}{2} \le Z \le \tfrac{1}{2}, \end{aligned}$$
(21)

and the deformed configuration is given by

$$\begin{aligned} r = r(R,Z), \qquad \theta = \Theta , \qquad z = z(R,Z), \end{aligned}$$
(22)

wherein \(R_{\text {in}}= \tfrac{R^*_{\text {in}}}{R^*_{\text {out}}}\) denotes the dimensionless inner radius. Of use in the sequel will be the aspect ratio of the undeformed airway, defined by \(\varepsilon = \tfrac{h^*}{R^*_{\text {out}}}\).

Under the above definitions, the dimensionless strain-energy functions for each constituent, and the whole tissue are given by

$$\begin{aligned} \Psi _c&= \frac{(I_1 - 3)}{2} + \frac{\alpha }{2} I_4, \end{aligned}$$
(23a)
$$\begin{aligned} \Psi _m&= \mu \frac{(I_1 - 3)}{2} + \frac{\omega }{2}\mathrm {H}(I_4 - 1)\left( \exp \left( \zeta (I_4 -1)^2\right) -1 \right) , \end{aligned}$$
(23b)
$$\begin{aligned} \Psi&= \Phi _c \Psi _c + \Phi _m \Psi _m, \end{aligned}$$
(23c)

where the dimensionless parameters \(\mu \), \(\omega \) and \(\alpha \) are defined by

$$\begin{aligned} \mu = \frac{{\mu ^*_m}}{{\mu ^*_c}}, \qquad \omega = \frac{\omega ^*}{{\mu ^*_c}}, \qquad {\alpha } = \frac{{\alpha }^*}{{\mu ^*_c}}. \end{aligned}$$
(24)

The dimensionless Cauchy stress tenors for each constituent, \(\pmb {\sigma }_i\), are then obtained from the dimensionless versions of (12); the dimensionless tissue stress is obtained from the dimensionless version of (15) or equivalently via \(\pmb {\sigma } = \sum _{i} \Phi _i{\pmb {\sigma }_i}\).

The dimensionless governing equations (6) and (16) (expressed in terms of the reference configuration) then read

$$\begin{aligned}&\frac{r}{R}\left( \frac{\partial r}{\partial R}\frac{\partial z}{\partial Z} -\frac{\partial r}{\partial Z}\frac{\partial z}{\partial R}\right) = 1, \end{aligned}$$
(25a)
$$\begin{aligned}&\frac{\partial r}{\partial R} \frac{\partial \sigma _{rz}}{\partial Z} - \frac{\partial r}{\partial Z} \frac{\partial \sigma _{rz}}{\partial R} + \varepsilon \left( \frac{\partial z}{\partial Z} \frac{\partial \sigma _{rr}}{\partial R} - \frac{\partial z}{\partial R} \frac{\partial \sigma _{rr}}{\partial Z}\right) \qquad \nonumber \\&\quad +\varepsilon \left( \frac{\partial r}{\partial R}\frac{\partial z}{\partial Z} - \frac{\partial r}{\partial Z}\frac{\partial z}{\partial R} \right) \left( \frac{\sigma _{rr} -\sigma _{\theta \theta }}{r}\right) = 0, \end{aligned}$$
(25b)
$$\begin{aligned}&\frac{\partial r}{\partial R} \frac{\partial \sigma _{zz}}{\partial Z} - \frac{\partial r}{\partial Z} \frac{\partial \sigma _{zz}}{\partial R} + \varepsilon \left( \frac{\partial z}{\partial Z} \frac{\partial \sigma _{rz}}{\partial R} - \frac{\partial z}{\partial R} \frac{\partial \sigma _{rz}}{\partial Z}\right) \qquad \nonumber \\&+\varepsilon \left( \frac{\partial r}{\partial R}\frac{\partial z}{\partial Z} - \frac{\partial r}{\partial Z}\frac{\partial z}{\partial R} \right) \frac{\sigma _{rz}}{r} = 0. \end{aligned}$$
(25c)

The dimensionless boundary conditions (17) and (18) are given by

$$\begin{aligned} r(1,Z)&= r_{\text {dis}}, \end{aligned}$$
(26a)
$$\begin{aligned} \pmb {\sigma }\left( R,\tfrac{1}{2}\right) {\mathbf {n}}_{\text {up}}&= {\mathbf {0}} ,\end{aligned}$$
(26b)
$$\begin{aligned} \pmb {\sigma }\left( R,-\tfrac{1}{2}\right) {\mathbf {n}}_{\text {low}}&= {\mathbf {0}} ,\end{aligned}$$
(26c)
$$\begin{aligned} \pmb {\sigma }\left( R_{\text {in}},Z\right) {\mathbf {n}}_{\text {in}}&= {\mathbf {0}} , \end{aligned}$$
(26d)

wherein \({\mathbf {n}}_{\mathrm{up}}\), \({\mathbf {n}}_{\mathrm{low}}\) and \({\mathbf {n}}_{\mathrm{in}}\) denote the the dimensionless unit normals to the upper, lower and inner surfaces, respectively.

Fig. 2
figure 2

Numerical results from the full thickness model (\(\varepsilon = 1\)) in various states of contraction (\(\alpha = 0\), \(\alpha = 0.1\) and \(\alpha = 0.2\)) and with a 5\(\%\) fixed stretch applied at the outer boundary of the PCLS. (i)–(iii) Radial deformation, r, and (iv)–(vi) axial deformation, z, plotted over the undeformed configuration, (R, Z). Cauchy stress components (vii)–(ix) \(\sigma _{rr}\), (x)–(xii) \(\sigma _{\theta \theta }\), (xiii)–(xv) \(\sigma _{zz}\) and (xvi)–(xviii) \(\sigma _{rz}\) plotted over the deformed configuration, (r, z). The parameter values are provided in Table 1 in Appendix C. Note that the colour bar scales differ between the individual plots

Full size image

3 Numerical results

Numerical solutions to (25)–(26) are obtained via the finite element method, implemented in the software FEBio (Maas et al. 2012). FEBio is a nonlinear finite element software specialising in computational biomechanics. In our application, we use FEBio to numerically solve the weak form of the conservation of linear momentum (i.e., Eq. (16)) assuming quasi-static equilibrium using linearised Newton-Raphson iterations. We simulated the deformation of a thin hollow cylinder (representing the PCLS) discretised by standard linear hexahedral elements (see Fig. 9 in Appendix A). Radial displacement matching experimental data was assigned to the outer radial surface (i.e., displacement boundary condition (17) applied with outward normal \({\mathbf {r}}\)), whilst the upper, lower and inner surfaces remain traction-free. We fixed z displacement for a single node on the outer boundary at \((R,Z) = (1,0)\) to eliminate the constant velocity solution. Figure 2 shows representative results, illustrating the mechanical response of the airway to an imposed radial stretch in the absence (passive case) and presence of active contraction. Details of convergence tests are given in Appendix A and parameter values, following that of Hill et al. (2018) or estimated from the literature, are given in the relevant figure captions (see also Table 1 in Appendix C). Although a consistent colour scheme has been used within all figures that follow, it should be noted that the scales differ between the individual plots (where indicated in figure captions) in order to capture the full qualitative behaviour of all of the results displayed.

In both the passive and active contraction cases we observe that the radial deformation, r, decays linearly with undeformed radius, R, but remains uniform axially (Fig. 2 (i)–(iii)). As required by the incompressibility of the material, the airway thins as it is stretched (Fig. 2 (iv)–(vi)).

The mechanical stress within the tissue displays significant spatial heterogeneity (e.g. Fig. 2 (xviii)). Moreover, we observe that while the deformation of the airway is qualitatively similar in the passive and active contraction cases (cf. Fig. 2 (i), (iii)), there are distinct qualitative and quantitative differences in the stress state between these regimes (cf. Fig. 2 (xvi), (xviii)). In particular, there is an increased and exaggerated heterogeneous stress distribution in the presence of active contraction. Furthermore, the axial dependence of these heterogeneous stress distributions increases with increased active contraction and is highlighted by the circumferential stress distribution, \(\sigma _{\theta \theta }\) (Fig. 2 (xii)).

In each case we observe increased radial stress, \(\sigma _{rr}\), at the outer boundary (in the direction of the prescribed stretch), with the stress at the inner wall remaining approximately zero in each case (Fig. 2 (vii)–(ix)). Similarly, tissue contractility significantly influences the circumferential stress, \(\sigma _{\theta \theta }\), as is to be expected, since the generated contractile stress acts in the direction of the circumferential fibres embedded in the airway (Fig. 2 (x)–(xii)). Moreover, we see that the circumferential stress is higher at the inner radius than at the outer radius (Fig. 2 (x)–(xii)).

Fig. 3
figure 3

The effect of reducing the PCLS aspect ratio, \(\varepsilon \), on the deformation and stress distributions across the radius of the airway wall. (i)–(ii) Radial deformation, r(R, 0), and (iii)–(iv) axial deformation, \(z(R,\tfrac{1}{2})\), plotted over the undeformed radius, R. Cauchy stress components (v)–(vi) \(\sigma _{rr}\), (vii)–(viii) \(\sigma _{\theta \theta }\), (ix)–(x) \(\sigma _{zz}\) and (xi)–(xii) \(\sigma _{rz}\) plotted over the deformed radius, r, at \(Z=0\). A 5\(\%\) fixed stretch is applied in the passive, \(\alpha = 0\), and active contraction case, \(\alpha = 0.2\). The aspect ratio decreases in direction of black arrows for \(\varepsilon \in \{1, 0.5, 0.25, 0.1, 0.05, 0.025, 0.01\}\) and the remaining parameter values are provided in Table 1 in Appendix C. Note that the scales of the y-axes differ between the individual plots

Full size image
Fig. 4
figure 4

The effect of reducing the PCLS aspect ratio, \(\varepsilon \), on the deformation and stress distributions through the axial thickness of the PCLS. (i)–(ii) Radial deformation, r, and (iii)–(iv) axial deformation, z, plotted over the undeformed thickness, Z. Cauchy stress components (v)–(vi) \(\sigma _{rr}\), (vii)–(viii) \(\sigma _{\theta \theta }\), (ix)–(x) \(\sigma _{zz}\) and (xi)–(xii) \(\sigma _{rz}\) plotted over the deformed thickness, z, at \(R=R_{\mathrm {mid}}\). A 5\(\%\) fixed stretch is applied in the passive, \(\alpha = 0\), and active contraction case, \(\alpha = 0.2\). The aspect ratio decreases in direction of black arrows for \(\varepsilon \in \{1, 0.5, 0.25, 0.1, 0.05, 0.025, 0.01\}\) and the remaining parameter values are provided in Table 1 in Appendix C. Note that the scales of the y-axes differ between the individual plots

Full size image

The thinning and stretching of the PCLS under the imposed stretch is reflected in the distributions of the axial, \(\sigma _{zz}\), and shear, \(\sigma _{rz}\), stresses, with an order of magnitude increase observed in the axial stress in the presence of contraction (Fig. 2 (xiii)–(xviii)). The positive (tensile) axial stress at the outer radius and negative (compressive) axial stress at the inner radius reflects the relative thickening at the inner radius compared with that at the outer. The shear stress is positive at the lower surface and negative at the upper surface reflecting the relatively increased displacement of material radially and downward at the upper (and upward at the lower) surface.

The preceding results correspond to an airway of thickness comparable to its outer radius (in particular, we set \(\varepsilon =1\)). The typical thickness for a PCLS is in the range of 100–\(250\mu \hbox {m}\) and a typical airway radius is in the range of 1–5mm, giving \(0.02< \varepsilon < 0.25\). Motivated by this, we investigate the dependence of the model behaviour on the PCLS thickness by varying the aspect ratio \(\varepsilon \). We consider the passive, \(\alpha =0\), and active contraction case, \(\alpha =0.2\), in Figs. 3 and 4 . Here, we reduce \(\varepsilon \) from \(\varepsilon = 1\) to \(\varepsilon = 0.01\) in the direction of the black arrows. Note that the variables are plotted as a function of deformed radius at the undeformed axial centre of the PCLS, i.e. at \(Z=0\) (Fig. 3). This is true in all cases apart from the axial deformation, z, which we plot as a function of radius at the undeformed upper surface of the PCLS, i.e. at \(Z=\tfrac{1}{2}\), in order to illustrate the thinning of the PCLS in response to stretch (Fig. 3 (iii), (iv)). The results illustrating axial dependence are all plotted over the thickness of the PCLS at the undeformed radial midpoint \(R = R_{\mathrm {mid}}\) (Fig. 4).

In both the passive and active contraction cases, the radial deformation at the axial centre line varies linearly with R and remains approximately invariant with \(\varepsilon \) (Fig. 3 (i), (ii); insets highlighting the very slight variation with \(\varepsilon \)) and with correspondingly little change in the radial stress (Fig. 3 (v), (vi)). Conversely, the near-uniform thinning of the airway, observed in Fig. 3 (iii), (iv), becomes marginally more exaggerated as \(\varepsilon \) is reduced and the PCLS thins more at the outer radius than at the inner radius. Similarly, the circumferential stress increases at the inner radius and decreases at the outer radius as \(\varepsilon \) reduces (Fig. 3 (vii), (viii)). On the other hand, we observe that the heterogeneous axial and shear stress distributions decay to zero in Fig. 3 (ix)–(xii).

The deformation and stress variation through the axial thickness is shown in Fig. 4. Here we observe only a weak dependence of the radial deformation and stresses on Z, that decays to uniformity with decreasing \(\varepsilon \). In particular, the axial stress decays to zero with \(\varepsilon \) (Fig. 4 (ix), (x)). In contrast, the axial deformation remains approximately linear in Z as \(\varepsilon \) decreases (Fig. 4 (iii), (iv); insets highlighting the very slight variation with \(\varepsilon \)). In the active contraction case, the above described features persist, but the variations in deformation and associated stresses are exaggerated quantitatively (cf. Figs. 3, 4).

4 Model reduction

Guided by the observations in Sect. 3, in this section we consider analytical simplifications of the biomechanical model (25)–(26). Firstly, in Sect. 4.1, we adopt a membrane model, following Wong and Shield (1969), that allows reduction to one spatial dimension. Subsequently in Sect. 4.2, we consider an asymptotic approach to obtain a reduced model describing the leading order PCLS deformation in the thin-PCLS-limit. In Sect. 4.3, we address the suitability of each approximation by comparing them to the full biomechanical model simulated in FEBio (Maas et al. 2012).

4.1 Membrane model

In this section we simplify the biomechanical model (25)–(26) by assuming that the PCLS behaves as an elastic membrane, in which case we neglect Z dependence and set \(z = Z\) so that there is no change in the axial thickness of the PCLS upon deformation. This membrane description has previously been used by Wong and Shield (1969) to approximate an axisymmetric stretch of an isotropic sheet; however, Wong and Shield (1969) find that the approximation breaks down when the sheet has an edge which is traction-free. The determinant of the deformation gradient approaches zero in the close vicinity of the traction-free edge and a singularity appears in the governing equations when the material is incompressible. Similarly, although we consider an anisotropic material with active contractile force generation, we find inconsistencies with the governing equations and the prescribed boundary conditions when the thickness of the PCLS is fixed. Specifically, the traction-free boundary conditions on the upper, lower and inner surfaces of the PCLS cannot be satisfied simultaneously, whilst preserving incompressibility, without the PCLS changing in thickness. Therefore, enforcing torsion-free axisymmetry as previously, we reduce the description of the PCLS to one spatial dimension and omit the traction-free boundary conditions on the upper and lower surfaces. The Cauchy stress, \(\pmb {\sigma }\), and pressure, \({\mathscr {P}}\), are functions of R only, satisfying

$$\begin{aligned} \frac{\mathrm {d} r}{\mathrm {d} R}&= \frac{R}{r}, \end{aligned}$$
(27a)
$$\begin{aligned} \frac{\mathrm {d} \sigma _{rr}}{\mathrm {d} R}&= \frac{2}{R}\left( {\psi }_1 + \psi _4 - \frac{{\psi }_1{R}^4}{{r}^4}\right) , \end{aligned}$$
(27b)

subject to the displacement outer boundary condition (26a),

$$\begin{aligned} r(1) = r_{\text {dis}}, \end{aligned}$$
(28)

and the free boundary condition at the inner radius (26d) (ommiting (26b)–(26c)),

$$\begin{aligned} \sigma _{rr}(R_{\text {in}}) = 0, \end{aligned}$$
(29)

and where the radial and circumferential stresses are constitutively defined as

$$\begin{aligned} \sigma _{rr}&= -{\mathscr {P}} + 2{\psi _1}\left( \frac{R}{r}\right) ^2, \end{aligned}$$
(30a)
$$\begin{aligned} \sigma _{\theta \theta }&= -{\mathscr {P}} + 2({\psi _1} + {\psi _4} )\left( \frac{r}{R}\right) ^2. \end{aligned}$$
(30b)

Integrating (27a) with respect to R, and imposing (28), we obtain

$$\begin{aligned} {r}^2 = {R}^2 - 1 + {r_{\text {dis}}}^2, \qquad {R_{\text {in}}} \le R \le 1. \end{aligned}$$
(31)

Subsequently, integrating (27b) with respect to R and applying the zero radial stress condition at the inner boundary (29) gives

$$\begin{aligned} \sigma _{rr}= \int _{R_{\text {in}}}^{R} \frac{2}{R'}\left( {\psi }_1 + \psi _4 - \frac{{\psi }_1{R'}^4}{{r}^4}\right) \mathrm {d}R'. \end{aligned}$$
(32)

In order to obtain \(\sigma _{\theta \theta }\) (30b), we require the pressure, \({\mathscr {P}}\); combining (30) and (32) provides

$$\begin{aligned} {\mathscr {P}} = 2\frac{{\psi }_1{R}^2}{r^2} - \int _{R_{\text {in}}}^{R} \frac{2}{R'}\left( {\psi }_1 + \psi _4 - \frac{{\psi }_1{R'}^4}{r^4}\right) \mathrm {d}R', \end{aligned}$$
(33)

and the constituent and total tissue stress follow directly.

4.2 Thin-PCLS-limit

In this section, motivated by the typical geometry of the PCLS (Sect. 3), we consider the limit \(0 < \varepsilon \ll 1\), so that the thickness of the PCLS is small in comparison to a typical airway radius. Correspondingly, and in view of our numerical results in Sect. 3 (in particular, Fig. 4 where we observe r becomes independent of Z for \(0 < \varepsilon \ll 1\)), we seek expansions of the form

$$\begin{aligned} r&= r^{(0)}(R) + \varepsilon r^{(1)}(R) + \varepsilon ^2 r^{(2)}(R) + {\mathcal {O}}(\varepsilon ^3), \end{aligned}$$
(34a)
$$\begin{aligned} z&= {z}^{(0)}(R,Z) + \varepsilon {z}^{(1)}(R,Z) + \varepsilon ^2 {z}^{(2)}(R,Z) + {\mathcal {O}}(\varepsilon ^3), \end{aligned}$$
(34b)
$$\begin{aligned} {\mathscr {P}}&= {\mathscr {P}}^{(0)}(R,Z) + \varepsilon {\mathscr {P}}^{(1)}(R,Z) + \varepsilon ^2 {\mathscr {P}}^{(2)}(R,Z) + {\mathcal {O}}(\varepsilon ^3), \end{aligned}$$
(34c)

adopting corresponding notation for the strain-energy functions where necessary and assuming that \(\Phi _c, \Phi _m, \omega , \zeta \) and \(\alpha \) all remain \({\mathcal {O}}(1)\) constants. We pause to highlight that the more general expansion, for which \(r = r(R,Z)\), and the leading term for \({\mathscr {P}}\) is \({\mathcal {O}}(\varepsilon ^2)\) (to obtain the proper leading order balance in the Cauchy stress), can be reduced to that shown in (34) (see Appendix C) and so we adopt this from the outset for brevity.

At leading order, the governing equations (25) read

$$\begin{aligned} \frac{\mathrm {d}r^{(0)}}{\mathrm {d} R}\frac{\partial z^{(0)}}{\partial Z}&= \frac{R}{r^{(0)}}, \end{aligned}$$
(35a)
$$\begin{aligned} 2 \frac{\mathrm {d}r^{(0)}}{\mathrm {d} R}\frac{\partial z^{(0)}}{\partial Z}\frac{\partial ^2 z^{(0)}}{\partial Z^2} - \frac{\mathrm {d} r^{(0)}}{\mathrm {d} R}\frac{\partial {\mathscr {P}}^{(0)}}{\partial Z}&= 0 , \end{aligned}$$
(35b)

subject to the displacement boundary condition at the outer radius (26a),

$$\begin{aligned} r^{(0)}(1) = r_{\text {dis}}, \end{aligned}$$
(36)

and the following free boundary conditions at the upper, lower and inner surfaces of the PCLS (26b)–(26d):

$$\begin{aligned} \frac{\mathrm {d} r^{(0)}}{\mathrm {d} R}\left( \left( \frac{\partial z^{(0)}}{\partial Z}\right) ^2 - \frac{{\mathscr {P}}^{(0)}}{2\psi _1}\right) = 0, \qquad \text {at~} Z = \pm \tfrac{1}{2}; \end{aligned}$$
(37a)
$$\begin{aligned} \frac{\partial z^{(0)}}{\partial Z}\left( \left( \frac{\mathrm {d} r^{(0)}}{\mathrm {d} R}\right) ^2 - \frac{{\mathscr {P}}^{(0)}}{2\psi _1}\right) = 0, \qquad \text {at~} R = R_{\text {in}}. \end{aligned}$$
(37b)

Equation (35a) and the boundary conditions (37a) provide

$$\begin{aligned} z^{(0)} = \lambda _z^{(0)}(R)Z, \end{aligned}$$
(38)

where the arbitrary function of R arising from the integration of (35a) vanishes due to axial symmetry in \(z^{(0)}\) about the axial centre line, \(Z=0\). Furthermore, the equation (37b) requires \({\mathscr {P}}^{(0)}={\mathscr {P}}^{(0)}(R)\). In view of which, together with the boundary conditions (37), we obtain:

$$\begin{aligned} {\mathscr {P}}^{(0)}(R)&= 2\psi _1{\lambda _z^{(0)}}^2(R), \end{aligned}$$
(39a)
$$\begin{aligned} {\lambda _z^{(0)}}(R_{\text {in}})&= \sqrt{ \frac{R_{\text {in}}}{r^{(0)}(R_{\text {in}})}}. \end{aligned}$$
(39b)

At \({\mathcal {O}}(\varepsilon )\) the linear momentum equations (25b)–(25c) read

$$\begin{aligned} \frac{\mathrm {d}r^{(0)}}{\mathrm {d}R}\left( \frac{\partial z^{(0)}}{\partial Z} \frac{\partial ^2 z^{(1)}}{\partial Z^2} - \frac{1}{4\psi _1} \frac{\partial {\mathscr {P}}^{(1)}}{\partial Z}\right) = 0, \end{aligned}$$
(40)

and the boundary conditions (26) provide

$$\begin{aligned} r^{(1)}(1)&= 0, \end{aligned}$$
(41a)
$$\begin{aligned} \frac{\mathrm {d}r^{(0)}}{\mathrm {d}R}\left( \frac{\partial z^{(0)}}{\partial Z} \frac{\partial z^{(1)}}{\partial Z} - \frac{{\mathscr {P}}^{(1)}}{4\psi _1}\right)&= 0, \qquad \text {at~} Z = \pm \tfrac{1}{2}, \end{aligned}$$
(41b)
$$\begin{aligned} \frac{\partial z^{(0)}}{\partial Z}\left( \frac{\mathrm {d}r^{(0)}}{\mathrm {d}R}\frac{\mathrm {d} r^{(1)}}{\mathrm {d} R} - \frac{{\mathscr {P}}^{(1)}}{4\psi _1} \right)&= 0, \qquad \text {at~} R = R_{\text {in}}. \end{aligned}$$
(41c)

Inspection of equation (41c) shows that \({\mathscr {P}}^{(1)}(R)\). In view of which, the boundary conditions (41b) and (41c) provide

$$\begin{aligned} {\mathscr {P}}^{(1)}(R)&= 4 \psi _1 \frac{\partial z^{(0)}}{\partial Z} \frac{\partial z^{(1)}}{\partial Z} , \end{aligned}$$
(42a)
$$\begin{aligned} {\mathscr {P}}^{(1)}(R_{\text {in}})&= 4\psi _1\left. \frac{\mathrm {d}r^{(0)}}{\mathrm {d}R}\right| _{R=R_{\text {in}}}\left. \frac{\mathrm {d} r^{(1)}}{\mathrm {d} R}\right| _{R=R_{\text {in}}}. \end{aligned}$$
(42b)

To summarise, we have reduced the problem to two leading order variables; the radial deformation, \(r^{(0)}(R)\), and the axial stretch, \(\lambda _z^{(0)}(R)\), and obtained the governing equation (35a) and the boundary conditions (36) and (39b). However, we require a second governing equation to determine the two variables, \(r^{(0)}(R)\) and \(\lambda _z^{(0)}(R)\). Therefore, we consider the \({\mathcal {O}}(\varepsilon ^2)\) momentum equations; (25b) reads

$$\begin{aligned}&\left( \frac{\partial ^2 z^{(0)}}{\partial R\partial Z} + \frac{R}{{r^{(0)}}^2}\right) \left( \frac{\mathrm {d}r^{(0)}}{\mathrm {d}R}\right) ^2 - \frac{1}{2\psi _1}\frac{\partial z^{(0)}}{\partial Z}\frac{\mathrm {d}{\mathscr {P}}^{(0)}}{\mathrm {d}R} - \frac{\psi _1 + \psi _4^{(0)}}{\psi _1R} + 2\left( \frac{\partial z^{(0)}}{\partial Z}\right) ^{-1}\nonumber \\&\quad \left( \frac{R}{{r^{(0)}}^2} -\frac{R^2}{{r^{(0)}}^3}\frac{\mathrm {d}r^{(0)}}{\mathrm {d}R} - \frac{R}{r^{(0)}}\frac{\mathrm {d}r^{(0)}}{\mathrm {d}R}\frac{\partial ^2 z^{(0)}}{\partial R\partial Z} \right) = 0. \end{aligned}$$
(43)

Equation (43) closes the leading order problem. Equation (25c) introduces higher order terms which are not of interest for the leading order problem and is therefore not needed here.

Substituting (38) and (39a) into (43) provides an equation for \(\lambda _z^{(0)}\) and hence, together with (35a), we obtain the following pair of coupled ODEs:

$$\begin{aligned}&\frac{\mathrm {d} r^{(0)}}{\mathrm {d} R} = \frac{R}{r^{(0)} {\lambda _z}^{(0)} }, \end{aligned}$$
(44a)
$$\begin{aligned}&\begin{aligned} \frac{\mathrm {d} {\lambda _z}^{(0)} }{\mathrm {d}R}&= \frac{R }{R^2 + 2{r^{(0)}}^2 {\lambda _z^{(0)}}^4}\left( 2{\lambda _z}^{(0)} - \frac{R^2}{{r^{(0)}}^2} \right. \\&\quad - \left. \left( 1 + \frac{\psi _4^{(0)}}{\psi _1}\right) \frac{{r^{(0)}}^2{\lambda _z^{(0)}}^2}{R^2} \right) . \end{aligned} \end{aligned}$$
(44b)

Together with the boundary condition (36) and the relation (39b), this provides a boundary value problem that describes the leading order radial and axial deformation. From this the leading order Cauchy stress components for the whole tissue and each of the constituents follow directly. We note that the boundary condition (39b) on \(\lambda _z^{(0)}\) is posed at the (unknown) deformed inner radius. We therefore solve (44) numerically by treating \(r^{(0)}(R_{\text {in}})\) as a shooting parameter and seek the solution set (\(r^{(0)}\), \(\lambda _z^{(0)}\)) at \(R_{\text {in}}\) that satisfies (36), (39b) and (44).

Fig. 5
figure 5

Comparison of the numerical simulations of the full model, with those of the membrane model, and the thin-PCLS-limit to demonstrate their validity. (i)–(ii) Radial displacement, r, (iii)–(iv) axial displacement, z, and (v)–(xii) Cauchy stress components, \(\pmb {\sigma }\), plotted as functions of undeformed radius, R, and deformed radius, r, at the undeformed centre of the PCLS (\(Z = 0\)), respectively. A 5% fixed stretch is applied to the PCLS in the passive, \(\alpha = 0\), (1st column) and active, \(\alpha = 0.2\), (2nd column) case. The aspect ratio \(\varepsilon = 0.01\) throughout and the remaining parameter values are provided in Table 1 in Appendix C. Note that the scales of the y-axes differ between the individual plots

Full size image

4.3 Suitability of approximations

In this section we compare numerical simulations of the full model (Eq. (25) with boundary conditions (26)), with those of the membrane model (Eq. (27) with boundary conditions (28)–(29)), and the thin-PCLS-limit (Eq. (44) with boundary conditions (36) and (39b)) to demonstrate their validity.

We plot the radial deformation, r, axial deformation, z, and the corresponding stresses, \(\pmb {\sigma }\), obtained in all three models, both in the presence and absence of active contractile force in Fig. 5. Results from the full and thin-PCLS models in Fig. 5 are plotted as functions of radius at the undeformed axial centre line (\(Z=0\)), apart from the axial deformation, z, which we plot as a function of radius at the undeformed upper surface of the PCLS (\(Z=\tfrac{1}{2}\)), in order to illustrate the thinning of the PCLS in response to stretch. Those from the membrane case, however, do not depend on Z; for illustrative purposes, we plot z fixed at \(Z=\tfrac{1}{2}\) in order to emphasise the thinning of the PCLS (relative to the reference configuration) that is displayed by the full and thin-PCLS models.

We find that, despite its relative mathematical and computational simplicity, the thin-PCLS-limit provides a suitable approximation to the full model, showing good quantitative agreement and excellent qualitative agreement in all variables. In contrast, the membrane model is unable to replicate the full model behaviour. The one-dimensional geometry of the membrane approximation constrains the inner radius of the membrane to deform corresponding to the displaced outer boundary in order to preserve incompressibility. As a result, we observe an increased radial deformation and elevated radial and circumferential stress in the membrane approximation compared to the thin-PCLS-limit approximation and the full model (Fig. 5 (i)–(ii), (v)–(viii), respectively). In contrast, the thickness of the PCLS in both the full model and the thin-PCLS-limit allows the generated stresses to be absorbed by the axial deformation (Fig. 5 (iii), (iv)). Hence, the thin-PCLS-limit provides a more realistic representation of the full problem (for small \(\varepsilon \)) than the membrane.

Active contraction accentuates the radial deformation of the PCLS and the thickness of the PCLS decreases accordingly in order to maintain tissue incompressibility (cf. Fig. 5 (iii), (iv)). This feature is only observed in the full model and the thin-PCLS-limit. Further contraction-induced deformation is not permitted in the membrane approximation due to the one-dimensional geometry and incompressibility constraint, and as a result, active contraction simply increases the stress generated in the membrane. Hence, there are significant qualitative and quantitative differences observed in the stress distributions of the two approximations and only the thin-PCLS-limit provides a suitable approximation to the full model.

Fig. 6
figure 6

The effect of stretch (applied at the outer boundary of the PCLS) on the constituent Cauchy stress components, \(\pmb {\sigma }_c\) and \(\pmb {\sigma }_m\), and the axial deformation at the upper surface, \(z^{(0)}(R,\tfrac{1}{2})\), as a function of deformed radius, \(r^{(0)}\), for different stiffnesses of ECM relative to that of ASM, \(\mu \). Simulation parameter values are provided in Table 2 in Appendix C. Note that the colour bar scales differ between the individual plots

Full size image

The full model exhibits rapid variation in the axial and shear stress components near the inner and outer airway radii (Fig. 5 (ix)–(xii)). However, this (small-amplitude) boundary layer behaviour in the axial and shear stress components near the airway boundaries (that is evident in the full model) is not captured by either of the simple models. Although the amplitude of these effects is very small, we believe that these are not numerical artefacts as they span multiple elements. Therefore, a boundary layer analysis of these features forms a natural future extension of this work.

5 Thin-PCLS-limit parameter exploration

In the preceding numerical experiments, our parameter choices follow that of Hill et al. (2018) or are estimated from the literature. The parameter values are provided in Table 2 in Appendix C. In this section we take advantage of the computational tractability of our reduced model ((36), (39b) and (44)) to explore the influence of the airway’s mechanical properties on the model behaviour, and in particular, examine differences in the constituents’ stresses that cause TGF-\(\beta \) activation. Such parameter exploration is computationally prohibitive in the full model.

The effect of the imposed radial displacement on the constituent stress, in the absence of contraction, is illustrated in Fig. 6. Here, we increase the fixed stretch applied from 0% (unstretched) to 20% (the maximum imposed in the experiments that are our primary motivation (Tatler 2016)). Over this range, we observe a significant increase in stress heterogeneity, with high radial (circumferential) stresses evident at the outer (inner) airway wall (Fig. 6). For both the ASM and ECM components, we see that the circumferential stress dominates over the radial stress (e.g., cf. Fig. 6 (i), (iv) and (vii), (x)). This is to be expected in the passive case due to the strain-stiffening of the ECM that occurs in direction of the circumferential fibres when stretched. However, we see that increasing the stiffness of ECM relative to that of ASM (\(\mu \)) dramatically alters the distribution of the constituent radial stress (cf. Fig. 6 (i), (iii) and (vii), (ix)).

As the imposed radial stretch is increased from zero, the axial deformation becomes less uniform across the airway thickness (Fig. 6 (xiii)–(xv)). Moreover, we observe that the heterogeneity of the axial thinning with increasing stretch is exaggerated with ECM compliance (cf. Fig. 6 (xiii), (xv)). This suggests that the stiffness of the ECM provides resistance to the imposed stretch across the airway wall, in addition to the associated strain-stiffening. As a result, the profile of the PCLS is more uniform for stiff ECM and thicker (thinner) at the inner (outer) radius for compliant ECM. When the ECM is relatively compliant, the stress state of the ASM is higher than that of the ECM for small stretches (due to the passive isotropic material properties, rm i.e. \(\mu \)) (cf. Fig. 6 (i), (vii)). However, strain-stiffening of the ECM increases exponentially with stretch and as a result, the stress state of the ECM increases more significantly with stretch than that of the ASM (cf. Fig. 6 (iii), (ix)).

Fig. 7
figure 7

The effect of constant contractile force, \(\alpha \), on the constituent Cauchy stress components, \(\pmb {\sigma }_c\) and \(\pmb {\sigma }_m\), and the axial deformation at the upper surface, \(z^{(0)}(R,\tfrac{1}{2})\), as a function of deformed radius, \(r^{(0)}\), in the presence and absence of fixed stretch applied at the outer boundary of the PCLS with amplitude 0%, 5% and 15%. Simulation parameter values are provided in Table 2 in Appendix C. Note that the colour bar scales differ between the individual plots

Full size image
Fig. 8
figure 8

The effect of the stiffness of ECM relative to that of ASM, \(\mu \), on the constituent Cauchy stress components, \(\pmb {\sigma }_c\) and \(\pmb {\sigma }_m\), and the axial deformation at the upper surface, \(z^{(0)}(R,\tfrac{1}{2})\), as a function of deformed radius, \(r^{(0)}\), in the presence and absence of a fixed stretch applied at the outer boundary of the PCLS with amplitude 0%, 5% and 15%. Simulation parameter values are provided in Table 2 in Appendix C. Note that the colour bar scales differ between the individual plots

Full size image

The influence of ASM contractility on the airway constituent stress and deformation, at fixed imposed stretch, is illustrated in Fig. 7. As expected, increasing the contractile force generated by the ASM leads to significant radial contraction of the airway, and associated resistance to axial thinning at the inner radius (Fig. 7 (xiii)–(xv)). Correspondingly, we observe elevated and increasingly non-uniform radial stress of the ASM and ECM constituents in Fig. 7 (i)–(iii) and (vii)–(ix), respectively.

In general, the stress distributions are qualitatively similar at each amplitude of fixed stretch, with a small stress increase at a larger fixed stretch. The circumferential stress of the ECM is an exception to this general observation and displays significantly different heterogeneous distributions for each imposed fixed stretch (cf. Fig. 7 (x), (xi) and (xii)). More specifically, in the absence of stretch, the circumferential stress of the ECM is maximal at the outer radius when the contractile force is high (Fig. 7 (x)). In contrast, in the presence of a 15% stretch, the circumferential stress of the ECM is maximal at the inner radius and when there is no contractile force (Fig. 7 (xii)). The transition between these two modes is evident in Fig. 7 (xi).

In contrast to our previous observations for increasing stretch in Fig. 6, we see that increasing the contractile force generated by the ASM leads to comparable radial and circumferential stress components of the ECM (Fig. 7). Here, the increasing contractile force and the strain-stiffened ECM leads to the constituents being comparably stressed.

The influence of constituent stiffness on the airway wall stress and deformation, at fixed imposed stretch, is illustrated in Fig. 8. The ratio of the passive isotropic stiffness of ECM to that of the ASM is given by \(\mu \). The ECM is more compliant than the ASM for \(\mu <1\), stiffer than the ASM for \(\mu >1\), and has the same stiffness as the ASM for \(\mu =1\). Here, we increase \(\mu \) in the presence of a constant contractile force generated by the ASM, \(\alpha \).

When the ECM is more compliant than the ASM (\(\mu <1\)) we observe a slight reduction in radial contraction compared to the case for which the ECM stiffness exceeds that of the ASM (\(\mu >1\)). Correspondingly, there is an increased resistance to axial thinning with increasing \(\mu \) observed in Fig. 8 (xiii)–(xv). As a result, we see that the magnitude of the stress components of the ASM decrease with increasing \(\mu \) (due to the decreased radial contraction), whilst the magnitude of the stress components of the ECM increase with increasing \(\mu \) (cf. Fig. 8 (v), (xi)). These observations persist and are emphasised under the application of fixed stretch (due to additional strain-stiffening of the ECM when stretched).

In general, we see that the stress distributions of the ASM display greater non-uniformity across the airway radius than that of the ECM, particularly when stretched. For example, the circumferential stress of the ASM and ECM differ significantly (cf. Fig. 8 (vi), (xii)). As the stiffness of the ECM increases, the circumferential stress of the ASM remains higher at the outer radius than at the inner radius in the absence of stretch (Fig. 8 (iv)). However, in the presence of a 15% stretch, the circumferential stress of the ASM varies significantly across the airway wall and is higher at the inner radius and lower at the outer radius (Fig. 8 (vi)). Therefore, imposed stretch induces a dramatic change in the distribution of the circumferential stress of the ASM and the apparent transition between these two extremes is observed in Fig. 8 (v). This behaviour is similar to that exhibited by the ECM when increasing the contractile force generated by the ASM in Fig. 7 (x)–(xii).

6 Discussion

Despite its prevalence in the population, the causes of asthma remain poorly understood; in particular, the feedback mechanisms linking inflammation, bronchoconstriction and cytokine activation are yet to be elucidated. To help address this, we develop a nonlinear fibre-reinforced biomechanical model of an airway in PCLS, an ex vivo assay widely used for studying asthmatic airway biomechanics. Our model accommodates agonist-induced ASM contractility and ECM strain-stiffening and allows us to examine the stress distributions of these individual constituents within the airway wall.

Direct numerical simulation of the model by means of the FEBio software (Maas et al. 2012) reveals the internal stress state of an axisymmetric airway within a PCLS under imposed deformation, and highlights the distinct qualitative and quantitative differences induced by ASM contraction. Such information is of key importance in interpreting PCLS experiments, and in particular those that seek to understand the above described feedback mechanisms. However, the computational complexity of this model precludes thorough investigation of the parameter space, mathematical analysis or coupling to time-dependence. To address this, we consider two reductions of the full model. First, we adopt a membrane representation, where axial variation is neglected a priori, and the PCLS thickness remains unchanged upon deformation. Secondly, and in view of the typical dimensions of the PCLS and numerical evidence, we consider an asymptotic reduction, appropriate for the limit in which the PCLS thickness is much smaller than the typical airway radius, and in which we are able to retain a description of the radial and axial airway deformation and the associated stresses. In each case, we reduce the model to one spatial dimension; the membrane model admits analytical solutions, while in the thin-PCLS-limit the model reduces to a pair of coupled nonlinear ODEs describing the deformation, numerical solutions to which are obtained via a shooting method. We find that the membrane model is unable to capture the full model behaviour, but that the asymptotically-reduced model provides a suitable approximation to the full model, at reduced computational cost.

Crucially, this computationally tractable model that we have developed allows for comprehensive investigation of the mechanisms underpinning pro-remodelling and contractile cytokine activation in asthmatic airways, a key aspect of the pathogenesis and presentation of asthma, that has only recently received attention. In particular, our future work will consider the positive mechanotransductive feedback loop between airway contraction and the activation of TGF-\(\beta \), that is implicated in long-term remodelling. Furthermore, we will include time-dependent loading of cyclic stretch to more faithfully represent the experimental protocol. Other important future considerations include developing our asymptotic reduction to accommodate, for example, spatially-dependent airway composition data (Brook et al. 2019) and undertaking parameterisation and model validation. In addition, as the PCLS thickness is reduced, we observe rapid variations in the axial and shear stress components near the inner and outer airway radii that are not captured by the asymptotic (or membrane) model; a boundary layer analysis of these features forms a natural extension of this work. More advanced theoretical considerations include development of a model accommodating unconstrained multiphase solid mechanics. This will allow for the differing strain rates of the two constituents to be examined, which may be important in understanding the cell-mediated activation of cyotkines in the PCLS.