A theoretical model of inflammation- and mechanotransduction-driven asthmatic airway remodelling

Inflammation, airway hyper-responsiveness and airway remodelling are well-established hallmarks of asthma, but their inter-relationships remain elusive. In order to obtain a better understanding of their inter-dependence, we develop a mechanochemical morphoelastic model of the airway wall accounting for local volume changes in airway smooth muscle (ASM) and extracellular matrix in response to transient inflammatory or contractile agonist challenges. We use constrained mixture theory, together with a multiplicative decomposition of growth from the elastic deformation, to model the airway wall as a nonlinear fibre-reinforced elastic cylinder. Local contractile agonist drives ASM cell contraction, generating mechanical stresses in the tissue that drive further release of mitogenic mediators and contractile agonists via underlying mechanotransductive signalling pathways. Our model predictions are consistent with previously described inflammation-induced remodelling within an axisymmetric airway geometry. Additionally, our simulations reveal novel mechanotransductive feedback by which hyper-responsive airways exhibit increased remodelling, for example, via stress-induced release of pro-mitogenic and pro-contractile cytokines. Simulation results also reveal emergence of a persistent contractile tone observed in asthmatics, via either a pathological mechanotransductive feedback loop, a failure to clear agonists from the tissue, or a combination of both. Furthermore, we identify various parameter combinations that may contribute to the existence of different asthma phenotypes, and we illustrate a combination of factors which may predispose severe asthmatics to fatal bronchospasms. Electronic supplementary material The online version of this article (10.1007/s10237-018-1037-4) contains supplementary material, which is available to authorized users.


Introduction
Asthma is a chronic lung disease characterized by inflammation, airway hyper-responsiveness (excessive bronchoconstriction in response to relatively low doses of contractile agonist; AHR) and airway remodelling. The last of these involves a series of structural changes, including thickening of the epithelial layer and subepithelial basement membrane (SBM; the collagen-dominated inner layer of the airway) and of airway smooth muscle (ASM) bundles (Holgate 2011;Brightling et al. 2012;James et al. 2012;Berair et al. 2013). While each of these three features contributes to asthma severity, how they interact is poorly understood. Most importantly, it is not clear whether AHR or remodelling are causes or consequences of the disease.
We hypothesize that while airway remodelling is initiated by inflammatory mediators, it is perpetuated by mechanical factors. The complexity of the underlying biochemical and mechanical processes, which span multiple length and timescales, makes identification of key interactions solely from biological experiments on isolated processes particularly challenging. Our aim is, therefore, to investigate the combined effect of repeated, short timescale, inflammatory episodes and associated mechanical forces, arising from ASM cell (ASMC) contraction, on long-term airway remodelling. To this end, we present a novel, quantitative mechanochemical modelling framework (informed by appropriate in vitro and in vivo studies) that integrates these processes for the first time. We use this model to elucidate emergent system dynamics and thereby identify key underlying pathogenic processes.
Although inflammation is considered to be the main process by which airway remodelling occurs, based on in vitro (Brightling et al. 2012;Dekkers et al. 2012;Noble et al. 2014 and references therein) and animal (Sjöberg et al. 2017;Alrifai et al. 2014;Silva et al. 2008;Zhu et al. 1999) studies, the causative effects of inflammation on remodelling are not fully supported by clinical trial or epidemiological data. For example, controlling inflammation with inhaled corticosteroids does not change the extent of airway remodelling or the decline of lung function with age (Guilbert et al. 2006;Strunk 2007). Moreover, airway remodelling may occur in early life in the absence of inflammation (James et al. 2012(James et al. , 2009. Bronchoconstriction, in the absence of inflammation, can also promote airway remodelling Oenema et al. 2013;Ge et al. 2012;Grainge et al. 2011;Tatler et al. 2011). In addition to the structural changes present in asthmatic airways, there is increasing evidence of altered baseline contractile tone that is thought to be the result of the chronic presence of contractile agonist or inflammatory mediators (Brightling et al. 2002).
It is not clear how this persistent tone arises, but it could enhance AHR (Bossé et al. 2009). Additionally, intra-subject heterogeneity of tone in asthmatics has been proposed as a cause of AHR (Brown and Togias 2016). Given that the area fraction of ASM in human biopsies has been shown to increase with the degree of asthma severity (Hassan et al. 2010), altered tone and ASM mass are likely interrelated. Furthermore, airway remodelling is also associated with extracellular matrix (ECM) changes (Kuo et al. 2011) and SBM thickening (Benayoun et al. 2003).
Inflammation in the airways is the protective response to allergen challenges, and is characterized by the recruitment of inflammatory cells (such as eosinophils), activation of resident mast cells, and over-expression of cytokines and chemokines (Brightling et al. 2003). Some of the latter interact directly with ASM to trigger contraction, and/or interact with mast cells causing them to degranulate, producing histamine and other contractile agonists (Kostenis and Ulven 2006). Inflammatory cells release mediators that also have the ability to induce remodelling (e.g. TGFβ; Halwani et al. 2011). Subsequently, the inflammatory cells and cytokines are gradually cleared from the tissue. Allergen challenges in humans are typically random, but in mouse models of asthma, inflammatory and contractile agonist challenges are administered and controlled artificially.
Tissue mechanics plays a significant role in airway remodelling and bronchoconstriction. For example, mechanical strain increases ASMC proliferation (Smith et al. 1994) and contractile protein expression (Smith et al. 1997(Smith et al. , 1995. In addition, TGF-β is a cytokine that mediates remodelling by inducing both cell proliferation and ECM protein production (Halwani et al. 2011), as well as having a potential contractile agonist role (Desmoulière et al. 2005;Grinnell and Ho 2002;Montesano and Orci 1988). It is activated by ASMCs (Tatler et al. 2011), likely during bronchoconstriction (Oenema et al. 2013;Tatler and Jenkins 2012) via mechanical stretch from latent complexes that are tethered to the ECM and to the ASMC (Froese et al. 2016;Noble et al. 2014;Wipff et al. 2007). Moreover, bronchial epithelial cells under mechanical compression shed growth factors such as endothelin-1 (ET-1), early growth response-1 (EGR-1) and TGF-β Drazen 2006, 2001;Tschumperlin et al. 2004) as well as to signal lung fibroblasts to express ECM proteins in vitro (Swartz et al. 2001).
The heterogeneous micro-mechanical stress environment, generated by ASM contraction, can be quantitatively and qualitatively very different in normal versus remodelled airways (Hiorns et al. 2016). Associated mechanotransductive processes may influence ECM deposition and ASM proliferation and migration. While it is recognized that these mechanisms play a crucial role in numerous developmental, physiological and pathological processes in other diseases (Hoffman et al. 2011;Martinez-Lemus et al. 2009), we have yet to understand how these mechanisms contribute to AHR and airway remodelling in asthma.
A significant body of work focusses on models of bronchoconstriction (e.g. Hiorns et al. 2014;Eskandari et al. 2013;Politi et al. 2010;Wang et al. 2008;Latourelle et al. 2002;Lambert and Paré 1997;Macklem 1996;Lambert et al. 1993), with much less attention paid to modelling of inflammation in the airways. Of particular relevance to this work are two of our previous models: firstly, our finitethickness continuum-based model ) of agonist-initiated contraction of the airway (represented as a nonlinear fibre-reinforced elastic material), which accounts for contractile force generation at the cell-level; secondly, our spatially-averaged (ordinary differential equation; ODE) model (Chernyavsky et al. 2014) of inflammation-driven switching of ASMCs from a contractile to a proliferative phenotype. The former predicts airway calibre changes and spatial stress heterogeneities, as a result of ASM contraction and in response to pressure fluctuations that mimic tidal breathing. In the latter, resolution of inflammation is implicated as the key factor in driving ASM mass accumulation; however, this model does not account for mechanotransductive feedback from mechanical stresses arising in the constricted remodelled airway wall. Thus, we combine and extend these two models by: (i) recasting our biomechanical model of the contractile airway into a well-established morphoelastic framework, that has been widely applied to growth and remodelling of soft tissues; and (ii) coupling this to the evolution of individual tissue constituents, using a multiphase description, governed by underlying biochemical processes such as inflammation-induced ASM phenotype switching, ASM cell proliferation or recruitment and ECM deposition, all of which can depend on mechanical stresses. Importantly, this coupled model accounts for two-way feedback between inflammation-driven changes in volume fractions of individual airway constituents and their resulting mechanochemical environment.

Mathematical formulation
We model the airway wall as a two-layered cylinder, representing the inner collagenous SBM and an outer (predominantly) smooth muscle layer, within which accumulation of ASM mass and ECM deposition is driven by both biochemical and mechanical stimuli (Figs. 1, 2), neglecting (for simplicity) the sub-mucosal epithelial layer. We first outline a description of the elastic deformation ( Fig. 1), consisting of the balance of linear momentum together with the constitutive specification of a hyperelastic mechanical law (Sect. 2.1.1). We then define the model inputs representing  Schematic of the evolution of the airway wall geometry with time, with the colours correspoding to tissue constituent volume fraction. Middle rows: Schematics of the airway geometry in the grown and stressed configurations, at the indicated time. At initial time t 0 , the airway wall is defined by inner radius R 1 , interface radius R int , and outer radius R 2 in κ 0 , the original, stress-free reference configuration. Airway wall growth is described by the mapping, G, from the original configuration to κ g , the grown, zero-stress configuration, defined by grown radii ξ 1 , ξ int (not shown), and ξ 2 at time t 1 . The airway is deformed via F to configuration κ, the current, stressed configuration, defined by radii r 1 , r int , and r 2 , subject to pressure boundary conditions P 1 and P 2 . Bottom row: Schematics of transmural mechanical stress distributions, generated from the active contraction and/or the elastic deformation, and contractile agonist. The mechanical stress, τ , modulates certain rates in the constitutive mass balance equations, e.g. (2.18, 2.21), to influence growth, G biochemical stimulation (Sect. 2.1.2). Finally (Sect. 2.1.3), we construct the mass balance equations, incorporating suitable biochemical and mechanotransductive processes (Fig. 2).

Geometry and kinematics
The airway tissue is modelled as a mixture (Ateshian 2007;Humphrey and Rajagopal 2002;Bowen 1976;Truesdell and Toupin 1960) consisting of four phases: ASMCs in either a contractile (c) or a proliferative ( p) phenotype; a collagendominated ECM (e); and an extracellular fluid (w) that also transports soluble nutrient (not modelled) for tissue maintenance. The outer layer is composed of multiple phases (predominantly ASMCs), and the inner layer, representing the SBM, is composed entirely of the ECM phase. Following the traditional continuum mechanics approach (Holzapfel 2000;Truesdell and Noll 1965), we assume initially (at time t 0 ) that each constituent a in the airway is in a (common) spatial, unstressed and unstrained reference configuration denoted κ 0 , in which the position of a particle is given by X.
We assume that the airway is an axisymmetric thick-walled cylinder of fixed length. We define fixed polar cylindrical co-ordinates R, Θ, Z in the radial, circumferential and axial directions, respectively, with the inner layer occupying R 1 R R int , and the outer layer R int R R 2 (Fig. 1). In Fig. 1, G is a topological mapping from κ 0 to an intermediate "grown" configuration κ g , in which the position of a particle originally at X is given by ξ (X, t), with cylindrical co-ordinates: where we have assumed the airway maintains axisymmetry and zero axial growth. The airway in the grown configuration κ g is deformed to a stressed configuration κ, and the position of a particle originally at ξ is now at x(ξ ), with the mapping given by the deformation gradient tensor F. The total deformation is given by using the standard multiplicative decomposition, H = FG. For concision, explicit dependence on time is suppressed here and throughout. For simplicity, we impose a plane-strain approximation, axisymmetry, and zero axial displacement. Thus, the deformation from κ g to κ is given by so that the elastic deformation gradient tensor is given by Assuming incompressibility, det(F) = 1 thus gives where r int = r (ξ int ). The initial airway geometry (Online Resource 2) is chosen to match the bovine airways used in LaPrad et al. (2010), from which the mechanical properties are obtained. These are similar in size and structure to generation 4 human airways (Harvey et al. 2013;Coxson et al. 2008;Williamson et al. 2011). The interface radius is chosen based on human airway histology (Benayoun et al. 2003). Mechanical properties We assume the tissue is a nonlinear hyperelastic heterogeneous (multiphase) anisotropic material (Bowen 1976;Truesdell and Noll 1965;Truesdell and Toupin 1960). The formulation for obtaining the constitutive mechanical relation for this type of material is given in detail by Ateshian and Ricken (2010) and Ateshian (2007). We neglect dissipative stresses and assume: that all solid constituents are constrained to move together (Humphrey and Rajagopal 2002); isothermality; electroneutrality; tissue Overview of the biochemical mechanisms. Allergen or contractile agonist challenges, f (t), specified to occur at times t i , drive evolution of an inflammatory factor, μ, and contractile agonist concentration, k. The magnitude and rate of clearance of μ and k are determined by constants a μ , a k , and c dμ , c dk , respectively. Inflammation leads to global release of contractile agonist at rate a kμ . Contractile agonist induces local ASMC contraction, and the resulting increased mechanical stress τ releases further cytokines, with contractile agonist properties, at rate a c . Contractile ASMCs undergo apoptosis (c a ) and switching to a proliferative phenotype (c cp ). The proliferative ASMCs divide (c p ) and switch to a contractile phenotype at a (high) constant rate (c pc ). Both inflammation and mechanical stress drive increases (from a baseline) in the contractile to proliferative switching rate (c cp (μ, τ )), and increasing mechanical stress drives increases in the proliferation or recruitment rate (c p (τ )). ECM proteins degrade (c de ) and are deposited (c be (μ)) at baseline rates during normal tissue maintenance, with the latter increasing with inflammation. Proliferating ASMCs produce additional ECM proteins (with rate c pe ). Blue dotted arrows indicate how constituent volume fractions are required for computation of the mechanical stress (τ ), as illustrated in Fig. 1. Rate constants are given in Online Resource 2 incompressibility; that stress arises from the elastic deformation only (i.e. the growth mapping does not impart stress); and that viscous stresses are negligible.
Following Hiorns et al. (2014), we apply the commonly used additive de-coupling of the active and passive Cauchy stress tensors: where 1 is the identity tensor and p is a Lagrange multiplier enforcing incompressibility, with the passive stress given by (2.6) in which we highlight that, since the mixture is constrained, F a = F. In (2.6), Ψ is the strain energy function per unit volume of the mixture, and W a is the strain energy function per unit constituent volume of constituent a (Ateshian 2007;Huyghe and Janssen 1997), noting that the fluid phase (w) does not contribute to the mechanical response of the tissue. The contractile ASMCs and collagen-dominant ECM form continuous fibre-like structures (Ijpma et al. 2017) and therefore are modelled as two sets of helical fibres wrapped symmetrically about the airway axis (to avoid torsion). The active stress in (2.5) is given by where A c is the contractile force density, defined as the force generated by the contractile ASMCs per unit area of constituent (Brook et al. 2010) and m ( j) c is the contractile ASMC fibre orientation vector (see Online Resource 1.1 for further details).
Specification of the constitutive mechanical response for the airway wall constituents In our previous work , the airway wall is modelled as a composite material of ASMCs and ECM, in which two families of fibres are embedded in an isotropic ground matrix, thus giving an anisotropic response. Fibre recruitment is modelled phenomenologically with an exponential dependence on stretch. Here, we de-couple the mechanical responses of the proliferating ASMCs, modelled as an isotropic neo-Hookean material; the contractile ASMCs, modelled as the fibreembedded material described above, and the ECM, modelled similarly under the additional assumption that collagen fibres bear load only after being extended beyond the recruitment stretch, λ u Hill et al. 2012;Robertson et al. 2011;Holzapfel and Ogden 2010). The strain energy functions for each of the constituents are given in Online Resource 1.2.
We introduce a scalar quantity τ representing the mechanical stress along the contractile ASMC fibre directions, given by which is used in the mass balance equations in Sect. 2.1.3 to elicit the mechanotransductive responses. The degree of contraction is directly related to the amount of agonist bound to the relevant contractile ASMC receptors. We therefore assume that the contractile force density, A c in (2.8), is a function of contractile agonist concentration, k, saturating as follows: where T c is a measure of hyper-responsiveness, K d represents the ratio of the dissociation rate of the ligand-receptor complex to its association rate, and n is the Hill coefficient, describing cooperativity. Balance of linear momentum In addition to the assumptions stated above, we further assume that body forces are negligible and inertial terms may be neglected due to slow timescales associated with quasi-static deformation. Therefore, in mechanical equilibrium, conservation of linear momentum requires ∇ · T = 0.
( 2.11) Under the assumption of no torsion and plane strain, (2.11) reduces to where T rr is the radial component and T θθ the circumferential component of the Cauchy stress. Boundary conditions for the elastic deformation For the elastic deformation, pressure boundary conditions are specified at the inner and outer radii, and continuity of radial displacements and stress at the interface, so that (2.13d)

Model inputs
We assume that a series of transitory allergen challenges drives a global inflammatory response that represent an influx of inflammatory cells, such as eosinophils, into the airway tissue, resulting in a cumulative inflammatory status denoted μ. The challenges may be administered (e.g. in chronic asthma mouse models), or occur naturally. Therefore, μ evolves according to where c dμ is the inflammation decay rate representing clearance of inflammatory cytokines or inflammatory cell apoptosis, and a μ the magnitude, and f (t) denotes the timing of events given by (2.16). Inflammation drives a local activation of mast cells and bronchoconstrictive mediators (represented by the concentration k) such as histamine or acetylcholine (Pelaia et al. 2008). The agonist induces active contraction of the ASM, which leads to airway narrowing together with associated local airway wall stresses. We assume that the local tensile stress, τ (2.9), induces activation of latent TGF-β; this cytokine also acts as a contractile agonist and therefore contributes to k. We also consider cases where local compressive stress drives remodelling. These cases represent compression-induced epithelial-cell-mediated expression of EGR-1 or ET-1, noting that, here, we do not model the epithelial cells directly.
As with inflammation, the frequency of contractile agonist challenges may be specified as model input to represent, for instance, artificial methacholine challenges in animals (Lauzon and Bates 2000; Gunst et al. 1988) or humans (Grainge et al. 2011). Contractile agonist, k, thus evolves over time according to where c dk is the agonist clearance rate, a k the magnitude of administered agonist stimuli, a kμ the rate of inflammationinduced agonist activation, and a c the rate at which the stress τ induces agonist release. The Heaviside step function H ensures that only tensile stresses release k (Wipff et al. 2007). By setting a μ = 0 or a k = 0 in 2.14 and 2.15, f (t) represents either contractile agonist-or inflammatory-only challenges, respectively. The challenges are represented by a series of Gaussian peaks where t i is a vector of N event times, and d and σ are constants (Chernyavsky et al. 2014).

Tissue growth
Under the assumption of axisymmetry, in cylindrical polar co-ordinates, the local density, ρ a , of each constituent (a = c, p, e, w) evolves according to: where ρ a and the (constrained mixture) radial growth velocity, v, are functions of the grown radius, ξ ( Fig. 1), and time, t; s a = s a ρ c , ρ p , ρ e represents the constituent-dependent source/sink terms, specified in detail below. We assume that the density of contractile ASMCs, ρ c , evolves through switching to and from a proliferative phenotype, ρ p ( Fig. 2; Naveed et al. 2017;Wright et al. 2013;Hirota et al. 2009). The definition of ASMC phenotype is based on the observable function of the cells arising from expression of intracellular proteins, e.g. proliferative ASMCs exhibiting decreased expression of contractile proteins (Wright et al. 2013). As in our previous model (Chernyavsky et al. 2014), the rate of switching is governed by the inflammatory status, μ, but here we additionally assume that switching can also be driven by the local fibre mechanical stress τ . Thus, where c pc and c a are positive constants, the first two terms representing switch back from the proliferative phenotype and apoptosis, respectively, the combination of which provides a logistic growth representation. c cp represents the inflammation-and stress-modulated rate of switching to the proliferative phenotype given by where c c0 , c c1 , c c2 and c cτ are positive constants. The Heaviside functions are used to divide the inflammation into three levels: healthy, mild, and severe, which are characterized by the thresholds μ 1 and μ 2 (Fig. 3). The final term represents local mechanical stress-induced switching: Fig. 3 Inflammation levels. Inflammatory status dynamics induced by a series of environmental stimuli, illustrating the parameters μ 1 , μ 2 , a μ , c dμ , and t i , noting that for periodic events, t i+1 − t i = Δt and ω = 1/Δt. The solid line represents relatively fast inflammatory clearance (high c dμ ), while the dotted line represents slow clearance Increases in proliferative ASMC density, ρ p , arise via proliferation and phenotype switching. Thus, proliferative ASMC turnover, s p , is represented by where the proliferation rate (c p (τ )) is modulated directly by the local stress τ according to ( 2.23) where the first term represents ECM synthesis by the proliferative ASMCs, e.g. via ASMC-mediated release of active TGF-β that induces ASMCs to synthesize collagen (Coutts et al. 2001); the second term represents modification of ECM by inflammatory mediators, e.g. via mast cell activation mediated by MMP-1; and the third term, linear degradation. The inflammation-driven ECM deposition is given by where, again, the Heaviside function is used to separate the three inflammation levels.
For simplicity, each of the constituents is considered intrinsically incompressible, i.e. their true densities (ρ T a ) remain constant in space and time (Ateshian 2011), and we assume that the true densities of the ASMCs and ECM are equal, a reasonable assumption in general (Gleason et al. 2004). In the following, we work in terms of the volume fraction, defined by Φ a = ρ a /ρ T a , and correspondingly write the source/sink terms as S a = s a /ρ T a . Assuming no voids, we obtain a Φ a = 1. (2.25) Equation (2.17), together with the definition of volume fraction and the source/sink terms, is re-expressed as where the volume fractions, Φ a , are functions of the grown radius, ξ (Fig. 1), and time, t, and a = c, p, e. As in Gleason et al. (2004), we assume that a constant and uniform tissue hydration is maintained, such that Φ w = 0.70. Summation of (2.26), together with (2.25), gives the mass balance equation for the entire mixture as We note that: Φ e is obtained using (2.25); S c , S p and S e are functions of Φ c , Φ p and Φ e . Thus (2.25), (2.27) and (2.26), together with initial conditions on Φ a (see Online Resource 2), ξ and boundary conditions on v, completely specify the time-evolving growth mapping G, provided that the mechanical stress state can be computed ( Fig. 1; Sect. 2.1.1).
Initial and boundary conditions on tissue growth During normal tissue maintenance, in the absence of inflammation or administered contractile agonist, we assume that ASM proliferation, recruitment and apoptosis, and ECM degradation and deposition in the airway wall balance to generate a homeostatic state. The non-trivial homogeneous steady state for (2.26), for which v = 0, is given by where we have made use of (2.18), (2.21), and (2.23). A linear stability analysis of (2.26) ensures a priori that the steady state (2.28) is stable. Hence, we impose this homeostatic steady state as the initial condition (2.29) A zero flux condition is imposed at the inner wall and at the interface between the two layers, which in cylindrical polar co-ordinates is given by In order to solve the ODE (2.27) representing growth in each layer, a boundary condition must be specified on the velocity. However, because (2.27) is first order, one is unable to specify the velocity at both the inner and outer boundaries of each layer. We therefore set the radial velocity to zero on the outer wall, so that all growth occurs inwardly (this choice is discussed further in Sect. 4). Additionally, we require the velocity and displacement at the interface of the two layers to be continuous. Hence, we have where the superscripts (i) and (o) denote limiting values taken from the inner and outer layers, respectively. Inflammatory or agonist challenge protocol Episodes of inflammation-inducing allergen or contractile agonist challenges are represented by (2.16) in (2.14) or (2.15), respectively. Simulations are performed over a 1000 day interval, with the challenges confined to the first 50 days, thereby allowing a long resolution period to investigate the effects of the challenges on long-term airway remodelling post-challenge. The numerical solution method is outlined in Online Resource 1.3.

Determination of material parameters
The passive material parameters are determined by fitting (2.12) to the quasi-static pressure-radius inflation data in LaPrad et al.

Results
To investigate the behaviour of our novel mechanochemical morphoelastic model, we first performed a one-at-a-time sensitivity study (Online Resource 3) on selected parameters of most relevance to airway remodelling. The analysis showed that remodelling was most sensitive to the magnitude and resolution rate of inflammation, during inflammatory challenges. However, due to the possible effect of combination of the large number of parameters in this model (see Online Resource 2) and the large differences in remodelling that occur between inflammation and contractile agonist challenges, we chose to perform a series of paired parameter explorations. These are used to investigate the effect of repeated inflammatory episodes, and mechanical forces that arise from repeated ASM contractions, on long-term airway remodelling and effective mechanical properties of the airway. Parameters are set to default values unless varied in the simulations. Results are discarded for parameter choices that lead to airway growth or contraction completely into the lumen, i.e. the inner radius decreases to zero. Unless otherwise specified, remodelling occurs in the outer layer only.
Effect of inflammatory challenges In our first set of simulations, we apply only inflammatory challenges (a k = 0), mimicking regular allergen exposures in experimental studies, and explore the response to changes in challenge frequency, ω μ , and inflammation resolution rate, c dμ . From an initial inner radius of 1.8mm, we observe inward airway remodelling towards the lumen, as depicted in Fig. 4a, showing the inner radius at 5 days post-final challenge for each parameter pair. In particular, we observe a "switch" effect, in which the response is insensitive to increases in ω μ and decreases in c dμ for a large region of parameter space, but dramatic increases in remodelling occur beyond a threshold parameter set. Contractile agonist retention time, defined as the number of days between the final inflammatory or agonist challenge and the reduction in the total amount of agonist in the airway cross-section to near zero (< 1x10 −6 ), remains rather low for relatively high ω μ and low c dμ , and is relatively insensitive to ω μ , but we observe a threshold in c dμ below which the agonist resolution time increases abruptly to 15 days from a baseline of approximately 10 days.
Detailed results for a specific parameter choice (highlighted by circles in Fig. 4a) are shown in Fig. 4b-e. Static pressure-radius curves are computed for time points during and following inflammatory challenges (Fig. 4b). For these parameter choices, the contractile agonist concentration remains low, since the rate of inflammation-induced agonist release is relatively low, and very little ASMC contraction occurs. During challenges, the pressure-radius curve shifts downward as the airway thickens due to inflammationinduced remodelling. Following challenges, the curves continue to shift slightly downwards, due to a small amount of post-challenge remodelling, prior to a return to the steady state.
The total amount of ASM and ECM increases during inflammatory challenges as ASMCs switch from contractile to proliferative phenotype, and their volume fractions relative to that of ECM increase over time (Fig. 4c, d). Relatively low mechanical stresses arise during challenges for this parameter set, as inflammation induces only a small amount of contractile agonist release, and hence contraction is also minimal. We observe compressive radial stresses in the airway mid-wall. In the circumferential direction, the stresses are found to be tensile towards the outer edge of the airway wall and compressive near the lumen (Fig. 4e). Mechanical stresses in the axial direction are compressive due to tissue incompressibility (Online Resource 5).
Effect of inflammatory challenge-induced ECM changes in the subepithelial basement membrane Next, to model SBM thickening associated with asthma, we allow for inflammation-induced ECM deposition and degradation in the SBM, as well as in the outer layer, by setting c e0 , c e1 , c e2 and c de in the inner layer to the default values given (for the outer layer). The addition of inflammation-induced ECM deposition results in increased inward remodelling (Fig. 5a),  Fig. 4, except here, inflammation also drives SBM thickening which, in turn, decreases effective airway compliance, as shown by modified pressure-radius curves (Figs. 5b,c,f,4b). Compressive circumferential stresses in the SBM are correspondingly reduced (Fig. 5c), due to increased cross-sectional area. In contrast, peak tensile circumferential stresses in the outer layers are greater with SBM thickening, presumably as a result of a thicker, stiffer inner layer and therefore an effectively stiffer airway.
Effect of contractile agonist challenges Here, we apply only contractile agonist challenges (a μ = 0), mimicking regular methacholine exposures in experimental animal studies. We highlight that, in these simulations, all remodelling is driven by mechanotransduction, i.e. by local, stress-mediated phenotype switching. We find that increasing the agonist challenge frequency (ω k ) and decreasing the contractile agonist resolution rate (c dk ) lead to increased remodelling at 5 days post-challenge. As in the inflammatory-only challenges, severe remodelling is observed only beyond a threshold parameter set. Contractile agonist resolution time increases with decreasing agonist resolution rate, c dk , with a threshold below which there is a dramatic change in resolution time from a baseline of approximately 5 days to over 60 days (Fig. 6a). The dramatic increase in resolution time is a result of the build-up of agonist concentration demonstrated by nonzero agonist concentrations at the end of each challenge (Online Resource 5). Similar to the inflammation-only challenges, the contractile agonist resolution time is relatively insensitive to the challenge frequency, ω k .
Despite the reduced ECM volume fraction in the outer part of the airway (see Online Resource 5), the overall effective stiffness of the remodelled airway is increased during contractile agonist challenges, as indicated by the continued downshift in pressure-radius curves at high transmural pressures (Fig. 6b). To inflate the airway, a relatively high transmural pressure is required to overcome the increased contractile forces in the strongly contracted state, with a strong downward shift at all pressures, and the appearance of a significantly more compliant portion that shifts towards the right. When the pressure is great enough to cause stretches that exceed the recruitment threshold (see Sect. 2.1.1 and Online Resource 1.2), there is a very rapid increase in the recruitment of ECM. Tissue growth (increased airway thick- ness) also contributes to the observed downward shift. The separation of the effects of contractile agonist and tissue growth becomes clear following challenges, where the compliant portion of the curve disappears (first curve following challenges, Fig. 6b), and the high pressure portion of the pressure-radius curves shifts downwards as remodelling continues post-challenge.
As contractile agonist challenges cause the airways to contract (Fig. 6c), local ASMC phenotype switching (Fig. 6d) and local release of contractile agonist increase with increasing tensile circumferential stress in the outer part of the airway wall (Fig. 6e) that arises from the agonist-induced contraction. Small amounts of locally activated residual contractile agonist continue to drive remodelling post-challenge, generating a feedback loop that slows agonist clearance (Fig. 6a), accounting for build-up in agonist concentration.

Effect of changes in intrinsic ASM hyper-responsiveness
To simulate the effects of increasing intrinsic responsiveness of ASMCs to contractile agonist, we return to the inflammation-only challenges and investigate the effect of paired changes in responsiveness (T c ) and challenge frequency (ω μ ). We find that remodelling increases with T c and ω μ (Fig. 7a), but a threshold effect only exists for increasing T c . Increasing T c also leads to much slower agonist clearance. At low values of T c , the agonist eventually clears (Fig. 7a, b), while at high values, a self-perpetuating feedback loop is established (Fig. 7a, c), due to the local mechanotransduction-driven agonist release, which is not resolved, and the airway eventually grows in to the lumen.
Global versus local effects of inflammatory and contractile agonist challenges In Figs. 4 and 6, we illustrate the effect of varying frequency and resolution rate of inflammation-and contractile agonist-only challenges, respectively. Here, we compare instead the effects of varying amplitude and resolution rate in these globally applied challenges (first and second columns, Fig. 8). Additionally, we compare these effects to changes in rates of inflammation-(a kμ ) and stress-mediated (a c ) local contractile agonist release (third column, Fig. 8).
Increasing the amplitude (a μ ) and decreasing the resolution rate (c dμ ) of inflammatory-only challenges lead to increased remodelling (Fig. 8a), without a clear threshold effect. Additionally, the degree of remodelling becomes less dependent on the inflammation amplitude for sufficiently fast resolution. At the (relatively low) default value of ASMC responsiveness, T c , employed here, the contractile agonist clears rapidly from the tissue upon cessation of inflammatory challenges, due to the small amounts of local agonist activated by mechanical stress (Fig. 8d). Low levels of agonist release during inflammatory challenges results in limited contraction and therefore very low mechanical stresses.
In contrast, we observe a very sharp threshold for increasing contractile agonist magnitude (a k ) and decreasing agonist clearance rate (c dk ) (Fig. 8b). Locally activated agonist remains in the tissue longer at lower clearance rates, and its effect is evident in the transmural variations in circumferential stress (Fig. 8f) in the time period following cessation of challenges (> day 50). High contractile agonist concentration induced by the agonist challenge generates significant bronchoconstriction and therefore relatively higher mechanical stresses than with inflammatory challenges.
Increasing the rate at which inflammation induces contractile agonist release (a kμ ) exacerbates remodelling, and increasing the rate of stress-induced agonist release (a c ) increases the positive mechanotransductive feedback, leading to additional remodelling (Fig. 8c). With both an increase in stress-mediated feedback and in inflammationinduced contractile agonist release, agonist resolution time Fig. 7 Effect of airway smooth muscle cell hyper-responsiveness on airway remodelling, active tone, and mechanotransductive feedback. Variation in a remodelled geometry with inflammation frequency (ω μ ) and hyper-responsiveness (T c ) parameter values. Transmural contractile agonist concentration is plotted as a function of radius at days 28, 30, and 32 for parameter value pairs indicated by the circled points on a: b where contractile agonist eventually clears from the tissue and c where contractile agonist remains in the tissue in an indefinite feedback loop, causing increasing remodelling long after cessation of challenges is increased, though it remains significantly faster than with direct contractile agonist challenge. Increasing a kμ and a c thus leads to an effective combination of inflammatory and contractile agonist challenges, as increased contrac- . 8 Global versus local effects on airway remodelling and active tone. Variation in remodelled geometry (1st row) with parameter pairs: inflammation magnitude and resolution (c dμ , a μ ; 1 st column), contractile agonist magnitude and resolution (c dk , a k ; 2 nd column), and inflammation-induced contractile agonist release and mechanotransductive agonist release (a kμ , a c ; 3rd column). Circumferential stress (2nd row) is plotted as a function of radius and time for parameter value pairs indicated by the circled points on the surfaces tion is observed during inflammatory challenges, and overall remodelling is higher than with agonist challenge alone. The higher levels of contractile agonist release during inflammatory challenges result in greater mechanical stresses (Figs. 8f) from agonist-induced bronchoconstriction. Effect of phenotype switching and proliferation rate modulation by tensile versus compressive mechanical stress In all of the above simulations, we have assumed that only tensile stresses can increase phenotype switching (via (2.19) and nonzero c f cp in (2.20)) and that the proliferation rate of the proliferative ASMCs is unaffected by local stress (default value c f p = 0). Here, we investigate the effect of varying these parameters (in inflammatory-only challenges) for both tensile and compressive stress-modulated phenotype switching and proliferation/recruitment rate. We vary either the stressinduced switching rate (c f cp ) or the proliferation rate (c f p ) and the contractile agonist clearance rate (c dk ). Transmural distributions of proliferating ASMC volume fraction at day 51 ( Fig. 9), selected from the overall results of this parameter exploration (Online Resource 6), illustrate the effects of these parameters on remodelling.
We observe similar amounts of remodelling for increases in both tensile and compressive stress-modulated pheno-type switching rate and decreasing agonist clearance rate with no clear threshold effect (see Online Resource 6). For the same parameter ranges, agonist resolution times are observed to be similar and relatively independent of c f cp for both cases. Distributions of the proliferative ASMC volume fraction significantly differ in the two cases, with larger volume fractions at the outer edge of the airway wall in the tensile stress-modulated case and at the inner edge in the compressive stress-modulated case (Fig. 9a), reflecting the observed circumferential stress heterogeneity and contractile agonist build-up. For our given initial conditions (Online Resource 2), both tensile and compressive stress-induced increase in phenotype switching results in a greater amount of airway remodelling than stress-induced increase in proliferation rate (Fig. 9a, cf. b). Agonist retention is similar between the two cases.

Discussion
The mechanisms underlying the interaction of inflammation, airway hyper-responsiveness and airway remodelling in asthma are poorly understood. Thus, we have devel- oped a new computational model of coupled biochemicaland mechanotransduction-induced remodelling of the airway wall. In the model, accumulation of ASM and ECM is driven by switching of ASMCs from a contractile to proliferative phenotype (which is capable of synthesizing ECM), and also, via a novel mechanotransductive feedback mechanism by which a mechanically stressed tissue releases mitogenic growth factors and contractile agonists.
Our results qualitatively match those reported for asthmatic airways (as well as results from our previous model; see Online Resource 7). For example, we have predicted narrowing of the lumen and a downward shift in the pressureradius curves that have been previously reported in humans (Williamson et al. 2011). Our model also predicts wall thickening of airways due to increased ASM volume, which has been reported both in asthmatic patients (James et al. 2012(James et al. , 2009) and in animal models (Alrifai et al. 2014).
The inflammatory twitch hypothesis suggests that selflimiting inflammatory events are invoked in the presence of an allergic stimulus and dissipate as the stimulus disappears, with long resolution periods being a possible explanation for the chronic inflammation characteristic of asthma (Pothen et al. 2016). Our results suggest that these long resolution periods may increase remodelling associated with severe asthma. In animal models of asthma, for example, transitory inflammatory cell recruitment, and increases in thickness of both ASM bundles and the SBM, have been associated with repeated/successive allergen challenges (Johnson et al. 2004;McMillan and Lloyd 2004). This response is thought to be the effect of pro-remodelling growth factors affecting activation of mast cells, ASM proliferation and ECM deposition (Naveed et al. 2017;James 2017) or increased recruitment of ASM progenitors such as myofibroblasts into the ASM bundles (Gerarduzzi and Battista 2017;Singh et al. 2014;Saunders et al. 2009). Thicknesses of both ASM bundles and the SBM have been shown to gradually return to control levels upon cessation of challenges (Alrifai et al. 2014;Leclere et al. 2012), with the resolution periods of inflammation and tissue growth being on the order of days to months (Pothen et al. 2016;Alrifai et al. 2014), which is reflected in our simulations, in which inflammatory status, μ, remains elevated with decreased clearance rate, c dμ .
Our work broadly follows that of Skalak (1980) and Rodriguez et al. (1994) in which growth and remodelling is treated as a topological mapping from a reference configuration to a grown configuration followed by an elastic deformation. A key component of this model is the implicit separation of timescales between the growth and the (instantaneous) elastic deformations. This approach has previously been applied to model airway narrowing (in which the folding that occurs during growth was modelled utilizing buckling theory; Moulton and Goriely 2011a, b), but this was in the context of a single phase, and growth kinematics were prescribed directly. Multiphase or mixture theory has been utilized in numerous cardiovascular studies (e.g. Valentin et al. 2013;Gleason and Humphrey 2005;Gleason et al. 2004;Rajagopal 2003, 2002), where the models were driven by evolution of the growth configuration of each phase, but without consideration of their interactions. To our knowledge, the approach taken here, in which the spatiotemporal evolution of the individual tissue constituents is driven by underlying biological and mechanochemical processes, has not been considered in airway remodelling. Thus, the multiphase model proposed here addresses the limitations of previous models by considering in detail the interactions between tissue constituents, as in the studies of Ateshian (2007) and Aparício et al. (2016). We use our model to carry out a series of parameter exploration studies to identify potential mechanisms underlying the pathogenesis and evolution of asthma, described below.
Impaired resolution of inflammation may explain airway remodelling characteristic of asthma As in our previous model (Chernyavsky et al. 2014), slower resolution of inflam-mation has a greater effect on airway wall thickening than either challenge frequency or amplitude (Figs. 4a, 8a), as a result of slower clearance of residual pro-remodelling and pro-contractile cytokines (See Online Resource 7). This effect may be responsible for increased remodelling and bronchoconstriction observed in patients with severe asthma that is poorly controlled with anti-inflammatory medications, e.g. corticosteroids. In these cases, inflammation may not be effectively cleared, e.g. by failure to induce inflammatory cell apoptosis (Wenzel 2012;Woolley et al. 1996).
Inflammation-independent bronchoconstriction-mediated mechanical stresses could drive airway remodelling Grainge et al. (2011) have demonstrated the possibility of inflammation-independent airway remodelling in asthmatics through application of methacholine-only challenges. Our model suggests the mechanisms by which this may occur. We show how tensile stresses along (circumferentially-oriented) muscle fibres, arising from agonist-induced contraction of ASMCs, can cause a local release of additional agonist (Fig. 6c). This mechanotransductive pathway may represent stretch-induced activation of latent TGF-β, which has both mitogenic (Halwani et al. 2011) and contractile agonist (Ojiaku et al. 2017;Oenema et al. 2013;Desmoulière et al. 2005;Grinnell and Ho 2002;Montesano and Orci 1988) properties. This additional local increase in agonist concentration triggers a local tensile stress-induced cell phenotype switching, thus driving increases in ASM mass. We have also shown that compressive stress-induced increases in phenotype switching or proliferation rate (Fig. 9), possibly mediated by shedding of growth factors such as EGFR and ET-1, could also explain the inflammation-independent remodelling observed in Grainge et al. (2011). Contractile agonists alone are insufficient to induce physiological changes in mice, as measured with plethysmography (Mailhot-Larouche et al. 2018). Our simulations suggest that very frequent challenges, impaired agonist clearance or increased rate of stress-driven agonist activation is required to cause significant remodelling (Fig. 8b). It is possible that an intrinsic inability to clear agonist, increased sensitivity of ASM/ECM to stress-driven TGF-β activation or EGFR/ET-1 shedding could place the asthmatic subject in a high-risk region of the parameter space.
Interaction of inflammation with intrinsic ASM hyperresponsiveness could explain persistent contractile tone and severe remodelling Our simulations show that increased intrinsic hyper-responsiveness causes increased remodelling (Fig. 7a) and significantly increased contractile agonist resolution times. This response is a result of increased bronchoconstriction at a given agonist concentration, driving increases in tensile stresses and hence further (mechanotransductive) activation of pro-mitogenic mediators and contractile agonists such as TGF-β. These results highlight that the retention of contractile agonist associated with increased hyper-responsiveness is a candidate mechanism accounting for the persistent tone observed in asthmatics (Brightling et al. 2002). As expected, long resolution times are associated with impaired clearance of contractile agonist, c dk (Fig. 6a). Notably, delayed contractile agonist resolution occurs when both the ASMCs are very hyper-responsive and inflammation challenge frequency is high. The combination of identified mechanisms could therefore be responsible for persistent contractile tone and remodelling, post-challenge, in hyper-responsive airways (Kariyawasam et al. 2007).
Increased airway remodelling could occur via local mechanotransductive effects An increase in both inflammationand mechanical stress-induced contractile agonist release (Fig. 8c) leads to a condition that is an effective combination of inflammation-only (Fig. 8a) and agonist-only challenge (Fig. 8b). Increased rates of inflammation-induced agonist release (a kμ ) could represent degranulation of the larger numbers of mast cells present in the ASM bundles of asthmatics (Naveed et al. 2017) and hence the production of contractile agonist (e.g. histamine). These released factors ultimately elicit a pro-mitogenic response where (tensile) mechanical stresses arising from ASMC contraction drive tissue remodelling via stress-induced phenotype switching of ASMCs. Separately, under the assumption of mechanically activated mediators (e.g. TGF-β), increasing stress-activated release (a c ) increases agonist-induced bronchoconstriction, leading to further agonist release and thus further contraction and stress-induced remodelling in a perpetual feedback loop. The feedback mechanism may explain how coupled inflammation and bronchoconstriction lead to increased remodelling in severe asthma.
Airway wall thickening may be a normal, local mechanoprotective response but ultimately becomes detrimental to global lung mechanics and function The increase in cross-sectional area from thickening of the SBM counterintuitively may serve to reduce mechanical stresses, arising during bronchoconstriction (Fig. 5c). This reduction in stress could, in turn, reduce compressive stress-driven shedding of growth factors mediated by epithelial cells. Our model does not allow for buckling of the SBM, due to our imposition of axisymmetry, but reduced compressive stresses as a result of the thicker SBM could also reduce propensity to buckle (Moulton and Goriely 2011a) and hence limit airway narrowing during bronchoconstriction.
Our results also indicate that remodelling-induced changes in functional mechanical properties of the airways may depend on what drove the remodelling. For instance, we show that the pressure-radius curves during inflammationonly challenges have qualitatively different characteristics to those associated with agonist-only challenges (cf. Fig. 4b, 6b).
Interaction of different underlying mechanisms may contribute to the existence of different asthma phenotypes A combination of various characteristics defines asthma subgroups, or phenotypes (see Wenzel 2012 and references therein). Here, we have identified various parameter combinations that may contribute to the existence of different asthma phenotypes. For example, the presence of eosinophils in the lung tissue is associated with thickened SBM, high expression of TGB-β, increased frequency and severity of symptoms, and more near-fatal events (Wenzel 2012;Miranda et al. 2004). Persistent eosinophils, despite corticosteriod treatment (typically associated with eosinophil apoptosis), is associated with adult-onset, less allergic severe (T H 2) asthma (Wenzel 2012). The slower clearance of inflammation leading to increased remodelling (Figs. 4a, 8a) and reduced contractile agonist clearance may underlie this type of asthma. Our results (Fig. 5) suggest that a thickened SBM may be associated with the increased presence of eosinophils, in agreement with the SBM thickness observed in non-eosinophilic compared with eosinophilic asthma (Miranda et al. 2004).
Our results show that for regular bronchoconstrictive events to induce long-term airway remodelling an extreme set of consequences (Figs. 6, 8b) is required, which may be the reason why, to our knowledge, contractile agonistinduced remodelling has not been observed in animal studies (e.g. Mailhot-Larouche et al. 2018). However, the study by Grainge et al. (2011) shows that contractile agonist challenges do lead to increased ASM thickness in humans. That mice do not spontaneously develop asthma but that humans do, suggests that the subjects of the Grainge et al. (2011) study may lie in a different region of parameter space, i.e. beyond the relevant thresholds.
Methacholine challenges have been shown to increase airway responsiveness in humans over the short term (Gazzola et al. 2017), and increased smooth muscle responsiveness may be independent of inflammation in some asthma phenotypes (e.g. non-T H 2 asthma; Wenzel 2012). Our model suggests a possible mechanism: with intrinsically hyperresponsive ASMCs, normal exposure to contractile agonists alone (in the absence of inflammation) may be sufficient to drive further agonist release (Fig. 6b), and increase ASMC responsiveness, thus leading an asthma symptom phenotypes not typically associated with inflammation, such as exercise-, obesity-induced asthma, or non-T H 2 asthma.
Increased ASMC hyper-responsiveness (Fig. 7) and/or increased mechanotransduction (Fig. 8c) may exacerbate remodelling, even at normal levels of inflammation (paucigranulocytic asthma; Carr et al. 2017). Thus, these mechanisms may underlie the degrees of severity associated with inflammation-induced (T H 2) asthma: mild to moderate asthmatics may be characterized by increased inflammation but normal ASMC responsiveness (Fig. 4), while severe asthmatics may have increased hyper-responsiveness that leads to pathological remodelling (Fig. 7). Alternatively, mechanisms that reverse influx of inflammatory cells, and hence, return inflammation to normal levels, may have become impaired (Pothen et al. 2016;Brightling et al. 2012). A remodelled airway exhibiting a combination of the characteristics described above may be primed for severe or even fatal bronchospasms in which airway contraction completely obstructs airflow.
Model limitations and future work Models of this nature necessarily require some simplifying assumptions. In particular, we have assumed that the airway remains axisymmetric (and free of shear) during both growth and elastic deformation. This simplified approach, in which only radial growth is considered, permits more straightforward analysis of the effects of varying concentrations of tissue constituents and contractile agonists across the airway wall. Nevertheless, our model is novel in accounting for local mechanotransductive effects and these spatial distributions. However, we note that some growth may occur circumferentially and, possibly, axially, in vivo. Circumferential growth/atrophy has been shown to underlie the development of residual stresses in arteries (Rodriguez et al. 1994); however, airways have been shown to exhibit little or no residual stresses (McKay et al. 2002), and so we omit this growth mode from the model for simplicity. Axial growth is similarly neglected, based on a lack of experimental data on axial growth of airways. The growth models of Ren (2013) and Grytsan et al. (2017) do account for both isotropic and varying degrees of anisotropic growth but do not consider heterogeneous spatial distributions. A more comprehensive approach would be to extend the models analysed herein to 3D. However, the numerical solution of the resulting equations is beyond the scope of the current study. Additionally, our imposition of axisymmetry prevents us from considering buckling of the stiff SBM. Future work could involve utilizing this model to predict regions of high compressive stresses in which buckling could occur (Moulton and Goriely 2011a).
To solve (2.27), we specify zero radial velocity at the outer wall, so that growth occurs inwardly. We note that this choice is somewhat arbitrary and that the growth velocity could be specified at any point along the wall thickness. We choose to specify zero growth on the outer boundary based on evidence that airways appear to narrow during asthma. Studies on human bronchial segments have shown that increased ASM in asthmatics contributes to exaggerated airway narrowing (Noble et al. 2013).
We do not allow for the separate evolution of the individual constituents' reference configurations, as in Humphrey and Rajagopal (2002). If our assumptions were relaxed so that these could move independently, interphase drag may be used to more accurately account for stress-induced activation of latent TGF-β. Moreover, changes in reference configurations, in addition to non-axisymmetric growth, could allow for the development of residual stresses. Although healthy airways exhibit negligible or small residual stresses (McKay et al. 2002), to our knowledge these have not been measured in remodelled airways. Relaxing this constraint would present a challenge in the theoretical development and in the numerical solution that is beyond the scope of the current work.
We have assumed that collagen fibre recruitment initiates at stretches above unity, with an exponential function representing gradually increasing fibre uncrimping. Rigorous experimental analysis of the airway structure-function relationship, through coupled mechanical testing and advanced microscopic imaging (e.g. Hill et al. 2012;Clifford et al. 2011), would enable incorporation of fibre dispersion and recruitment to refine our model (e.g. Gasser et al. 2006;Sacks 2003;Lanir 1983Lanir , 1979. Also, we have not accounted for diffusion of applied contractile agonists, or of activated cytokines such as TGF-β. Finally, migration of myofibroblasts (ASM progenitors) into the muscle bundle has been hypothesized as a mechanism of increased ASM. Although this migration (and differentiation into ASMCs) has not been modelled here explicitly, the ASMC proliferation in our model could instead represent this recruitment of ASM progenitors, and the timescale associated with migration is implicitly accounted for in the proliferation rate (reinterpreted as a 'recruitment' rate of myofibroblasts which then become ASM).
Notwithstanding these limitations, key advantages of our model are the ability to generate measurable biological and mechanical output that may be tested experimentally and the ability to separate the long-term effects of growth from the relatively shorter term mechanical effects of pressurization and active contraction. A study, combining the modelling approach described here with experimental measurements of tissue constitution, geometry, and active mechanics, is required to distinguish between the effects of geometric remodelling and changes in tissue mechanics.
Therefore, this work forms the theoretical basis of a model that is currently being tested against data from in vivo animal experiments, utilizing an ovalbumin model of inflammationinduced asthma in mice. We have shown how the volume fraction of constituents may be evaluated at any time along the remodelling process (e.g. Fig. 4). These will be compared to measurements made on histological sections taken from animal models of asthma, to infer underlying mechanisms as discussed above, and also for further model development and validation. Additionally, we have demonstrated the ability to assess mechanical properties of the airway as they evolve, via computed pressure-radius curves (e.g. Fig. 4b) which are commonly measured experimentally, to show changes in passive and active tissue mechanics and distinguish between healthy and diseased airways.
In addition to airway remodelling, the approach outlined in this paper has broad applicability to other areas of tissue mechanics. For example, the mechanotransductive feedback mechanisms described herein likely have direct applications to models of aneurysm growth (Gryt-san et al. 2017;Aparício et al. 2016). Moreover, the coupled biochemical-biomechanical elements are relevant to inflammation-associated adventitial collagen deposition observed in arteries under hypertension (Bersi et al. 2016). Furthermore, the technique we present for modelling a non-homogeneous spatial distribution of tissue constituents would be applicable to the myocardium, since local changes in tissue constitution occur during mechanical overload, e.g. in hypertension-induced myocyte hypertrophy and collagen deposition (Hill et al. 2014) and during remodelling following myocardial infarction (Gajarsa and Kloner 2011).
Conclusions Our results suggest that mechanical stresses, arising from bronchoconstriction, initiated by multiple inflammatory or contractile agonist challenges, and driving agonist release, generate a mechanotransductive feedback loop. With this feedback loop, increased ASM hyperresponsiveness to contractile agonists and impaired clearance of inflammatory factors and contractile agonists leads to increased remodelling, increased bronchoconstriction, and maintenance of an increased baseline contractile tone due to the chronic presence of contractile agonists. The key factors that allow such a state to emerge are (i) increased intrinsic hyper-responsiveness of ASMCs to contractile agonists, (ii) delayed contractile agonist clearance from the tissue, (iii) increased release of contractile agonist by inflammatory cells, and (iv) mechanotransductive feedback, in which active contraction generates increases in mechanical stress that subsequently initiate the release of additional contractile agonists and thus exacerbate remodelling. Targeting mechanotransductive pathways (e.g. TGF-β) and rapid resolution of inflammation or contractile agonist may have potential for treating severe forms of asthma.