Coupled Immunological and Biomechanical Model of Emphysema Progression

Chronic Obstructive Pulmonary Disease (COPD) is a disabling respiratory pathology, with a high prevalence and a significant economic and social cost. It is characterized by different clinical phenotypes with different risk profiles. Detecting the correct phenotype, especially for the emphysema subtype, and predicting the risk of major exacerbations are key elements in order to deliver more effective treatments. However, emphysema onset and progression are influenced by a complex interaction between the immune system and the mechanical properties of biological tissue. The former causes chronic inflammation and tissue remodeling. The latter influences the effective resistance or appropriate mechanical response of the lung tissue to repeated breathing cycles. In this work we present a multi-scale model of both aspects, coupling Finite Element (FE) and Agent Based (AB) techniques that we would like to use to predict the onset and progression of emphysema in patients. The AB part is based on existing biological models of inflammation and immunological response as a set of coupled non-linear differential equations. The FE part simulates the biomechanical effects of repeated strain on the biological tissue. We devise a strategy to couple the discrete biological model at the molecular /cellular level and the biomechanical finite element simulations at the tissue level. We tested our implementation on a public emphysema image database and found that it can indeed simulate the evolution of clinical image biomarkers during disease progression.

Introduction

Chronic obstructive pulmonary disease (COPD) is estimated to affect more than 500 million people worldwide, causing significant disability, loss of quality of life and social burden, with costs in excess of € 56 billion per year in the European Union (Decramer et al., 2012). The disease has a lifetime prevalence of about 28% and cigarette smoking is commonly considered to be the principal risk factor (Gershon et al., 2011). Recent projections suggest that COPD will be the third cause of global mortality by the year 2030.

The pathogenesis of COPD is still not completely understood (Larsson, 2007; Yoshida and Tuder, 2007) and involves a number of multi-scale cellular processes, including airways inflammation, adaptation and innate immunity to cigarette smoking, sensitivity to self and not-self antigens, accelerated senescence, and deregulation of mechanisms of cell repair (Repapi, 2010; Pavord et al., 2012). Interactions between the environment and a selected group of candidate genes is also considered very important (Akinbami et al., 2012; Mizuno et al., 2017; Zhao et al., 2017).

Clinical management of COPD involves consistent use of inhaled corticosteroids that help reducing COPD mortality. However, their efficacy is limited (Faner and Agustí, 2016) and many patients experience exacerbations and poor symptoms control (Brightling et al., 2012).

As a matter of fact, the clinic presentation of COPD is not homogeneous, but presents two main clinical phenotypes, emphysema and chronic bronquitis, each with many sub-types, different comorbidities and risk profiles (Martinez et al., 2012). Even if there is no therapeutic target that can reverse the decline of lung function over time (Vestbo et al., 2013), a broader recognition of markers associated with adverse risk (Partridge et al., 2006) and therapies that specifically target different phenotypes specifically reduce exacerbations and improve patient's life (Castro et al., 2010; Holgate, 2012).

An additional problem is that it is extremely challenging to do an early detection and staging of COPD. This is because the gold standard for clinical diagnosis is Pulmonary Function Tests (PFTs) which is not sensitive enough to detect any disease progression before a large part of the lung has been compromised (Cooper et al., 2017). It is also not sensitive enough to detect different subtypes and elucidate different mechanisms of actions.

Specifically in emphysema, the continued inflammation of lung parenchyma eventually leads to a loss of collagen and elastin in the alveoli (Sharafkhaneh et al., 2008; Goldklang and Stockley, 2016). As a result of this sustained damage the septa become increasingly compliant and eventually fail mechanically during normal breathing. This reduces the area available for gas exchange causing dyspnea and shortness of breath. In addition, the mechanical damage due to emphysema is likely to stimulate tissue repair mechanisms at cellular level, that result in the production of type I collagen (Crosby and Waters, 2010). As a matter of fact, alveolar fibrosis is observed in emphysematous spaces, in the form of thickened and stiffened alveoli, which most likely contributes to shortness of breath (Yousem, 2006).

Faced by this complexity in the mechanisms and the lack of a simple clinical tests, it is important to assess the patient by integrating information from heterogeneous sources such as molecular data and medical imaging, in order to adapt the treatment options with the phenotype and risk profile. A promising option is to include information from computational models of biological systems that can account for causative effects, otherwise difficult to apprehend in clinics. These models have the potential to predict complex behaviors, elucidate regulatory mechanisms, and inform experimental designs to eventually point out specific factors to control or therapeutic targets, in order to improve patient management (Di Ventura et al., 2006).

Cancer research has already exploited computational models over different spatial and temporal scales as a promising way to describe complex diseases (Deisboeck et al., 2009, 2011; Wang et al., 2015). There, multiscale models interact with clinical data to generate and test different hypotheses, facilitating drug development (Clancy et al., 2016) and optimizing delivery and therapeutic effect (Cristini et al., 2017). We refer the interested reader to the detailed review by Wang and Maini (2017).

Recent interdisciplinary advances contributed to unravel the complex pathophysiological mechanisms that occur in COPD on both the macroscopic and microscopic scale. In case of macroscopic model of the respiratory system, for example Bordas et al. (2015) describes how to obtain a specific mesh of the patient for CFD simulations and Berger et al. (2016) discuss the application of a poroelastic deformation model for pulmonary ventilation. Chernyavsky et al. (2014) proposes a theoretical model of the possible effect of inflammation on the restriction of small airways. The reader can also refer to the review of COPD multi-scale modeling by Burrowes et al. (2013).

Among others, the “Protective Artificial Respiration” initiative fundamentally contributed to the understanding of COPD. We would like to cite Wiechert et al. (2011) for their multiscale model of respiratory system that coupled large bronchi and small alveoli, as well as Roth et al. (2017a,b) respectively for a study of the essential interactions between flow and deformation in the lungs and a simplified model of lung microstructures. Also Verdugo et al. (2017) reported on efficient solvers for respiratory mechanics. Among the works devoted to particle deposition we recall the work of Freitas and Schröder (2008) for a numerical study of 3D flows in a human lung model, and Lintermann and Schröder (2017) for the simulation of aerosol particle deposition and Calmet et al. (2016) for their model and simulations of particle deposition based on High-Performance Computing. Very recently, an experimental characterization of the nonlinear compressible behavior of the parenchyma is reported in Birzle et al. (2018).

For the microscopic modeling, literature contains numerous works on the modeling of the immune system at the molecular level. For instance, Folcik et al. (2007) developed an agent-based model for the innate and adaptive immune system while (An, 2008) contributed an agent-based model of the epithelium. A model of inflammation with interactions between macrophages and fibroblasts capable of simulating scarring, tissue damage and fibrosis is presented in Brown et al. (2011). Most of the studies on AB modeling of COPD focus on emphysema, and mainly study the resulting destruction of the tissue. The most common method uses a 2D network of springs to represent alveolar tissue (Mishima et al., 1999). These modeling studies have the merit to highlight the redistribution of forces within the tissue during the progression of emphysema. This simulated progression was found to produce experimentally observed emphysema patterns (Suki et al., 2003) and was extended to 3D by Parameswaran et al. (2011) through the use of cuboidal cells to represent the alveoli. The European AirProm project has initiated the study of multi-scale models for the study of COPD (Burrowes et al., 2013).

In the INSPIRE project¹, we would like to give a multi-scale, multi-physics description of the phenomena that cause the onset of emphysema and the possibility to predict the risk profile of the patient. Accordingly, the main purpose of the presented work is to propose a multi-scale model, able to integrate known interactions among inflammation, remodeling and parenchyma destruction, with particular attention to the role played by the immune system. We extended our previous work Ceresa et al. (2017) to couple the dynamics of the biological events captured through agent cooperation in an agent-based (AB) model with a biomechanical simulation of the tissue captured by a coupled Finite Element (FE) Model that iteratively predicts the evolution of the mechanical cues transmitted to the cells inside the lungs. We hope that such model could, once properly refined and validated, add to the interpretation of the specific disease phenotype toward the prediction of personalized risk profiles. We think our model builds nicely on the previous cited literature for the microscopic models because we use a less simplified model of the molecular interactions. In addition, we explicitly take into account the mechanical forces the tissue is subjected to using a well-vetted FE model, while others have worked more on connection models with elastic spring.

In the following sections we will discuss the coupled AB and FE model that we contribute (section Methods), and the experimental setup designed to validate the model on a public CT dataset of emphysema images (section Experimental Setup). We then present the results of the experiments, their discussion (section Results and Discussion) and the conclusions and future works (section Conclusions and Future Works).

Methods

As we commented before, research and clinical practice suggest that emphysema development happens along two different time-scales: a slow molecular one due to the inflammatory response to solid particles (Cosio et al., 2009), and a rapid one, caused by sudden rupture of the alveolar walls due to mechanical forces which act on lung tissue during respiration (Suki et al., 2003).

In the following sections, first we present a dynamic model of inflammatory response using ordinary differential equations (ODE) taken from literature that does not account for spatial and mechanical effects (section Well-Mixed Molecular Model of Inflammation and Tissue Remodeling). This is followed by an AB molecular model for inflammation and remodeling coupled with a FE model of biomechanical tissue that supersedes those limitations (section Agent Based Model of Inflammation and Coupling to the Finite Element Model).

Well-Mixed Molecular Model of Inflammation and Tissue Remodeling

In order to prepare the implementation of the AB model and define the rules thereof, we performed a large bibliographical study to obtain relevant information about:

• cytokines IL₁, IL₈, IL₁₀, TNFα and TGFβ production,

• macrophage migration, activation and differentiation into M₁ and M₂ types,

• feedback loops in the production of pro- and anti-inflammatory cytokines

• the role of MMPs on collagen cleavage and fibroblast deposition which are important terms for elastin degradation and remodeling.

This literature (Ignotz and Massagué, 1986; Onozaki et al., 1988; Oliver et al., 1993; Bellingan et al., 1996; Tsutsumi et al., 1996; Darby et al., 1997; Meng and Lowell, 1997; Hehenberger et al., 1998; Horio et al., 1998; Cobbold and Sherratt, 2000; Steinmüller et al., 2000; Eberhardt et al., 2002; Huang et al., 2002; Maass et al., 2002; Zhang et al., 2003; Mantovani et al., 2004; Porcheray et al., 2005; Tanaka et al., 2005; Edwards et al., 2006; Lenga et al., 2008; Marino et al., 2008; Moro et al., 2008; Jin and Lindsey, 2010; Wang et al., 2012) is reported in Reference section and it is associated to the different biological parameteres considered in Table 1. We focus mainly on the well-vetted interactions between different types of macrophages, pro- and anti-inflammatory cytokines, fibroblasts, collagen deposition and degradation, neutrophils and elastase production. Those interactions were already described by Brown et al. (2011), Jin et al. (2011), and Wang et al. (2012) and our main contribution was to integrate all the available information of the different biological processes and adapt them for the specific case of emphysema modeling. The final model we used is composed by two algebraic equations and thirteen coupled non-linear ordinary differential equations (ODE). This model belongs to the category of well-mixed (WM) systems in the sense that no spatial effects are considered.

TABLE 1

Table 1. Parameters of the AB model.

These equations are presented below (Equations 1–15) and the biology they reflect can be schematically represented in an integrated picture of the main molecular and cellular actors that regulate the chronic immune response and the consequent changes in tissue properties (Figure 1), after initial particle deposition on the lung tissue.

FIGURE 1

Figure 1. Agent based model of tissue destruction in emphysema progression. Particles coming from inhaled smoke cause secretion of cytokines such as TNFα and TGFβ by the epithelial cells. Those act, at first, as chemotactic factors and attract undifferentiated alveolar macrophages and fibroblasts. and the alveolar macrophages. Secondly, they induce the activation of the macrophages and their differentiation in the M₁ and M₂ subtypes. Those will create a delicate dynamical balance between inflammatory and anti-inflammatory signals such as IL₁, IL₈, and IL₁₀ that affect the activation of protease such as MMPs, the recruitment of neutrophils and fibroblasts. MMPs directly cleavage the collagen from the tissue and are responsible for the deposition of abnormal collage that leads to fibrosis, together with fibroblasts. Elastase destroys the elastin in the tissue. Both abnormal collagen deposition and reduction in elastin deteriorate the mechanical properties of the tissue.

The aforementioned particle deposition causes sustained inflammation of the tissues with a fast secretion of Tumor Necrosis Factor alpha (TFNα-Tα in the equations for brevity) and a slow secretion of Transforming Growth Factor beta (TGFβ-Tβ in the equations for brevity) respectively by monocytes (M) and epithelial cells. These cytokines attract monocytes, according to the model proposed by Wahl et al. (1987) (Equation 1):

$\begin{array}{ll}M(T_{\alpha })=0.335T_{\alpha }^3-6.309T_{\alpha }^2+32.281T_{\alpha }+57.302& (1)\end{array}$

and further govern the differentiation between inactivated macrophages (M_un) and the specific sub-types M₁ and M₂, according to Equations (2–4):

$\begin{array}{l}{\dot{M}}_{un}=M(T_{\alpha })-k_2{M}_{un}\frac{I{L}_1}{I{L}_1+c_{IL1}}-k_3{M}_{un}\frac{T_{\alpha }}{T_{\alpha }+c_{T\alpha }}\\ -k_4{M}_{un}\frac{I{L}_{10}}{I{L}_{10}+c_{IL10}}-\mu M_{un}& (2)\end{array}$

$\begin{array}{l}{\dot{M}}_1=k_2{M}_{un}\frac{I{L}_1}{I{L}_1+c_{IL1}}+k_3{M}_{un}\frac{T_{\alpha }}{T_{\alpha }+c_{T\alpha }}+k_{m21}{M}_2\\ -k_{m12}{M}_1-\mu M_1& (3)\end{array}$

$\begin{array}{ll}{\dot{M}}_2=k_4{M}_{un}\frac{I{L}_{10}}{I{L}_{10}+c_{IL10}}-k_{m21}{M}_2+k_{m12}{M}_1-\mu M_2& (4)\end{array}$

We see that all attracted monocytes will become inactivated macrophages first, and then switch to one of the two sub-types depending on constants k₂₋₄ and the concentration of pro-inflammatory cytokines IL₁, TFNα and anti-inflammatory cytokine IL₁₀. The pro-inflammatory cytokines will promote differentiation to M1 and the anti-inflammatory ones to M2. The promotion effect of the cytokines is mediated by the Hill equation for a cooperative binding type (Stefan and Le Novère, 2013) with coefficients c_IL1, c_Tα, and c_IL10. In time, the transition from M₁ to M₂ can be reversed, with constants k_m12 and k_m21. Eventually, the macrophages will be removed through the lymphatic system with rate μ.

Continuing our discussion of Figure 1, we see that each macrophage type will now secrete cytokines with a dynamic expressed by Equations (5–7):

$\begin{array}{ll}{\stackrel{.}{IL}}_{10}=k_5{M}_2\frac{c_1}{I{L}_{10}+c_1}-d_{IL10}I{L}_{10}& (5)\end{array}$

$\begin{array}{ll}{\dot{T}}_{\alpha }=k_6{M}_1\frac{c_1}{I{L}_{10}+c_1}-d_{T\alpha }{T}_{\alpha }& (6)\end{array}$

$\begin{array}{ll}{\stackrel{.}{IL}}_1=k_7{M}_1\frac{c_1}{I{L}_{10}+c_1}-d_{IL1}I{L}_1& (7)\end{array}$

Here we see in Equation (5) that IL₁₀ is secreted by M₂ proportionally to k₅ and regulated by self-inhibition with effectiveness c1. Eventually, it is degraded with half-time decay rate d_IL10.

In Equations (6, 7) we have an analogous process for the secretion of TFNα and IL₁ by the M1 macrophages subtype.

Additional TGFβ is secreted from fibroblasts (F) and M₂ to increment deposition of collagen in the composition:

$\begin{array}{ll}{\dot{T}}_{\beta }=k_8{M}_2+k_{9}F-d_{T\beta }{T}_{\beta }& (8)\end{array}$

$\begin{array}{ll}{F}_g(T_{\beta })=0.05T_{\beta }^3-0.98T_{\beta }^2+6.54T_{\beta }+7.11& (9)\end{array}$

$\begin{array}{ll}\dot{F}=k_{10}{F}_g(T_{\beta })F-d_{F}F& (10)\end{array}$

$\begin{array}{ll}\dot{C}=k_{11}F-d_{FC}MMPC& (11)\end{array}$

$\begin{array}{ll}\stackrel{.}{MMP}=k_{12}{M}_1-d_{FC}MMPC-d_{M}MMP& (12)\end{array}$

where K₈ is the secretion rate by M₂, k₉ the one by fibroblasts and d_Tb the decay rate in Equation (8). Fibroblasts proliferates from the population of already existing cells proportionally to TGF_β in Equations (9–10) and emigrate with rate d_f. Collagen deposition is governed by Equation (11), where we have to consider the deposition rate k_FC, and the degradation effect of matrix-metallo-proteinases (MMP). Those are enzymes produced by M₁ that degrade the collagen, as described in Equation (12).

Finally, macrophages attract neutrophils to the wound site by secreting IL₈, and those release the elastase enzyme that cleaves the elastin bonds in the fibers:

$\begin{array}{ll}{\stackrel{.}{IL}}_8=k_{13}{M}_2\frac{c_2}{I{L}_8+c_2}-d_{IL8}I{L}_8& (13)\end{array}$

$\begin{array}{ll}\dot{N}=k_{14}(1-\frac{N}{N_{\max }})\frac{I{L}_8}{I{L}_8+c_{IL8}}-{\mu }_{N}N& (14)\end{array}$

$\begin{array}{ll}\dot{E}=k_{15}N-d_{E}E& (15)\end{array}$

IL8 secretion (Equation 13) is similar to Equation (5), with constant k₁₃, a self-inhibition term with efficacy c₂ and a degradation constant of d_IL8. Equation (14) governs the recruitment of neutrophils up until their maximum value N_max with a cooperative effect of IL₈ and an emigration rate of μ_N. Finally the density of elastase is dependent upon the number of neutrophils and the inactivation rate, d_E.

The final proportion of elastase and collagen density is directly used in our biomechanical model to calculate the properties of the lung tissue for the FEM simulation as discussed at the end of the next section.

All the values for the discussed parameters are presented in Table 1 and the related literature is listed in Reference section.

Agent Based Model of Inflammation and Coupling to the Finite Element Model

In order to add spatial effects to the molecular model of inflammation and tissue remodeling, an AB model is created, using Equations (1–15) as a basis for the behavior of the agents. The first important difference is that the simulation of the agents happens on a grid. This gives the model an inherent spatial aspect and allow us to consider additional details w.r.t. the WB model. For instance, now the composition of the alveolar unit (AU) -which includes among others epithelial cells, collagen, elastin and basement membrane (Zemans et al., 2015)- becomes relevant. In our case, every cell of the grid represents a small portion of the AU with different variables accounting for the content of elastin, collagen, the cytokines and the structural integrity of the cells (called “tissue-life” in the following).

During the simulation a “smoking” signal determines whether we introduce particles into the simulated AU or not. This signal is a periodic square wave with frequency f _s and intensity e_s. The intensity quantifies the exposure, that is, the number of particles inhaled in each cycle. The signal starts from zero and last for a total smoking time of T_s. By varying frequency, intensity and total time, we can study the effect of particles on the model as detailed in experiment of section Experiment to Characterize Parameter Sensitivity. After the end of the total smoking time, the model is allowed to run for some additional time steps in order to reach equilibrium again.

The initial, unperturbed, dynamic of the system includes a small number of inactivated macrophages that move randomly, “patrolling” the tissue and searching for solid particles, similarly to the mononuclear cells behavior described by Auffray et al. (2007). When the smoking signal is active, inhaled particles deposit and cause an initial rapid rise of TFNα that attracts inactive macrophages to the deposition site according to Equation (2). From there, according to the dynamics described in Equations (3–4), macrophages differentiate in M1 or M2 subtypes which respectively govern the production of pro-inflammatory (Equations 6, 7, 12) and anti-inflammatory cytokines (Equations 5, 8, 13).

As previously indicated in Brown et al. (2011), at sites with high levels of pro-inflammatory cytokines, tissue is damaged by a complex network of interconnected factors called Damage-associated Molecular Patterns (DAMPS) (Matzinger, 2002; Lotze et al., 2007). This aspect was not included in the WM model because of its specific spatial nature, but it is implemented in the AB model where the tissue life of the AU is reduced proportionally to the inflammation level. Damaged tissue (i.e., with reduced tissue-life) in turn, start secreting TFGβ to recruit fibroblasts for wound healing as in Equations (9, 10).

The model tracks separately the amount of collagen and elastin in the tissue and their equilibrium varies depending on the concentration of fibroblasts, neutrophils, elastase and MMPs as in Equations (11, 12, 14, 15).

The cellular death caused by DAMPS and the amount of collagen and elastin, all affect the mechanical properties of the tissue used in the FE simulations. In this first version we use an elastic, isotropic material implemented in Elmer FEM software (Råback, 2013). Now, on the one hand, when a cell dies, we reduce its Young's modulus (E_Tissue) to 1 Pa, to account for the fact that it contributes no more to the elastic properties, but without changing the topology of the mesh. On the other hand, if the cell is not dead, its Young's modulus is calculated as a linear mixture of the corresponding concentration of elastin and collagen as in Equation (16). The initial values are E_el = 0.1 kPa and E_cl = 20 kPa, as described by Suki et al. (2011).

$\begin{array}{l}{E}_{Tissue}={\delta }_{el}{c}_{el}{E}_{el}+{\delta }_{cl}{c}_{cl}{E}_{cl}\\ {\delta }_{el}=0.7;{\delta }_{cl}=0.3,c_{el}\in [0,1];c_{cl}\in [0,1];& (16)\end{array}$

The material properties are calculated and loaded in the solver as continuous static field using a custom made Fortran code.

Apart from the molecular damages caused by DAMPS, the tissue can also die because it was subjected to too much strain during the mechanical simulations. While elastin withstands deformation as high as 100%, the maximum tensile strain of pure collagen fibrils with low cross-link density is considered to be around 10% of the initial length (Depalle et al., 2015; Sherman et al., 2015). Accordingly, the maximum tensile strain for each cell is calculated weighting the previous values for the amount of elastin and collagen contained.

We present in Figure 2 the indirect coupling strategy used for the AB and FE models. At each step the former simulate additional particle deposition that accounts for continued smoking; release of inflammatory cytokines and degradation of mechanical properties. Periodically the AB model is frozen and the calculated tissue properties are imported in the AB-FE coupler code which will reconstruct a topologically equivalent geometry, recover the contours of the damaged zones and assign new material properties taking into account the final amount of collagen and elastin from the AB model. The resulting information is passed to the FE solver that runs until convergence and then export the strain results for further processing. After the FE solver has run, the second coupler code, FE-AB is run to import the strain field and calculate which fibers, if any, have been destroyed in the simulation. It thus updates the AB status and restarts it with the updated state.

FIGURE 2

Figure 2. Full coupled model. Original patches from a public emphysema database are segmented to separate the parenchyma from the vessels and airways and seed deposition is simulated. For each seeded pixel and for all its neighbors we run a simulation job that represent the evolution of 130 alveoli. In each job there is a cyclic sequence between the agent and finite element model. At each step the former simulate additional particle deposition that accounts for continued smoking; release of inflammatory cytokines and degradation of mechanical properties. Periodically the AB model is frozen and the calculated tissue properties are imported in the AB-FE coupler code which will reconstruct a topologically equivalent geometry, recover the contours of the damaged zones and assign new material properties taking into account the final amount of collagen and elastin from the AB model. The resulting information is passed to the FE solver that runs until convergence and then export the strain results for further processing. After the FE solver has run, the second coupler code, FE-AB is run to import the strain field and calculate which fibers, if any, have been destroyed in the simulation. It thus updates the AB status and restarts it with the updated state.

Experimental Setup

The next sections deal with the more experimental part of our work. First, we explain the inner working of the coupling between AB and FEM solvers (section Procedure to Couple AB and FEM). Then, we present the meshing process (section Mesh Creation and Sensitivity). In the central part of this section we detail the two main experiments that validate our implementation: the first is an initial exploration of the sensitivity of the model to initial parameters (section Experiment to Characterize Parameter Sensitivity), while the second is the validation on a public CT image dataset of emphysematous lungs (section Experiment to Study the Emphysema Progression in Clinical Images). Finally, we briefly discuss the High Performance Computing infrastructure we used to run the studies (section High Performance Computing).

Procedure to Couple AB and FEM

In a typical execution cycle, the AB model is stopped at regular intervals and control is transferred to the FE model for analysis of the mechanical strains. After each interruption of the AB simulation, the latest iteration of this simulation is saved to disk and the AB-FE coupling code first calculates the percentage of damaged tissue area, as predicted by the AB model, and evaluates whether there is enough healthy tissue to proceed with the mechanic simulation. If this is the case, the saved status of the AB model is inspected to retrieve the last topology of the computational grid and the amount of collagen and elastin is used to calculate the new Young modules of the tissue according to Equation (16). This information is used together with connected component analysis, morphological operators and k-Nearest Neighbors (kNN) classifiers, to extract the contours of the broken tissue, define a 2D mesh and assign mechanical properties to each element. Materials, boundary conditions and solver parameters are adjusted if necessary and a case directory is created for the FE solver. The FE model runs asynchronously until convergence of the steady state and deformation and displacement fields are saved in a vtk compatible format (vtu). After that, the FE-AB coupling code is executed again. It reads back the strain fields from the solver status files and determines which, if any, nodes of the mesh have exceeded their maximum strain. Those are added to the damaged zone and the agent simulation is restarted. Cycle by cycle the coupled simulations continue until tissue damage is above 80% of the area or until the desired simulated time is reached.

A detailed view of typical results for the inflammation, meshing and mechanical process is shown in Figures 3, 4.

FIGURE 3

Figure 3. Progression of parenchyma destruction in Agent Based Model. From left to right we see the effect of increased inflammation, tissue damage and final destruction. (A) Concentration of proteases in a small sample of the tissue during model execution. (B) Due to the continued effect of the high proteases levels the tissue is damaged. (C) Snapshot of the damaged tissue as sent to the FEM model. Tissue in foreground (white) has greatly diminished mechanical properties.

FIGURE 4

Figure 4. Detail of the meshing process for a patch with initial emphysema formation. (A) The geometry is automatically generated and meshed by gmsh using our code taking into account Agent Based Model. (B) The resulting mesh has an adaptive size to ensure fast convergence and is able to capture complex shape for the destroyed tissue. (C) The generated mesh is then connected to our FEM solver and simulated until convergence is reached. In this image, the displacement field is shown.

Mesh Creation and Sensitivity

We use 2D FE meshes with topologies equivalent to the AB simulation grids. The exact size and topology of each mesh is thus very dependent on the current state of the simulation. In addition, an optimization step is run after the first mesh creation using Gmsh². The average mesh contains around 50,000 polygons with four nodes. We manually refined the parameters of the mesh creation to ensure a quick convergence of the FE simulations, while keeping a low computational cost, necessary to ensure reasonably fast and smooth interactions between the AB and the FE models.

Experiment to Characterize Parameter Sensitivity

We studied several possible parameters to characterize the model's behavior. First, we varied the quantity of particles and the frequency with which they are added. We varied the number of particles inhaled in each smoking step from 0 to 20 and the smoking time from 10 to 90 simulation steps, for a total of 25 experiments. This will be referred to as the “exposure” experiment in the results section.

In order to asses the sensitivity to the parameters, we selected and varied six main parameters of the model as described in Table 2. This resulted in a total of 2⁶⁻² = 16 experiments following a fractional factorial analysis. We will refer to this as the “parameters” experiment in the results section.

TABLE 2

Table 2. Parameters experiments.

Experiment to Study the Emphysema Progression in Clinical Images

One of the main objective of this model was to predict the development of emphysema in time. We devised an initial way to test our hypothesis using a public lung image dataset. We explain our approach in the following sections.

Dataset

We test our system against the public CT Emphysema database (Sorensen et al., 2010). We use 168 square patches manually annotated in a subset of the 115 high-resolution CT (HRCT) slices. As explained in the previous reference, CT scanning was performed using General Electric (GE) equipment (LightSpeed QX/i; GE Medical Systems, Milwaukee, WI, USA) with four detector rows. The acquisition protocol was: in-plane resolution 0.78 × 0.78 mm, slice thickness 1.25 mm, tube voltage 140 kV, and tube current 200 mAs. The slices were reconstructed by using a high-spatial-resolution (bone) algorithm. The data comes from a study group of 39 subjects, including 9 never-smokers, 10 smokers, and 20 smokers with COPD. Figure 2 shows a sample of each of the three categories of images.

Pre-processing

All slices were automatically segmented and reviewed to create a mask of only parenchyma tissue. In order to prepare the computational model, we first segmented the pulmonary tissue in the lung patches, using a fixed threshold of −750 HU. Stereological analysis of the lung parenchyma revealed a mean of 500 million alveoli per double lung in the normal population, with a mean alveolar volume of around 4.2 × 10⁶ μm³ and, on average, 170 alveoli per cubic millimeter (Ochs et al., 2004). In our case, with an anisotropic spacing of 0.78 × 0.78 × 1.25 mm³, this corresponds to roughly 130 alveoli per voxel. For each voxel of this binary mask we generate a planar grid of the 130 alveoli that is used as a computational mesh.

Particle Deposition and Simulation

As detailed in Figure 2, once the patch has been segmented, random pixels of the parenchyma and their neighbors are marked as “affected” and, for each one, a new simulation of the AB and FE models is run. Final results are mapped back into the main image patch and the updated mechanical properties calculated by the coupled AB-FE model are linearly translated back into HU values. In this way, we simulate the typical darkening of the CT scan caused by emphysema progression.

High Performance Computing

Among many different frameworks available for AB modeling (Abar et al., 2017), we chose to use Pandora (Rubio-Campillo, 2014), for its ease of programming and superb scalability. The model is implemented in an in-house version, specifically modified to allow biological model developments and available online³.

In order to satisfy the high demand in computational resources, we run the simulations on our institution's supercomputing SNOW Linux cluster. The cluster is currently composed by 20 computing nodes and a total of 840 cores with a theoretical calculation capacity of 8.49 Tflops. Highly relevant for agents simulations were six GeForce GTX TITAN X GPU with 12 Gb of memory.

Results and Discussion

In this section we present the results of the two main experiments that we have used to validate our implementation. Those experiments were previously discussed in detail respectively in sections Experiment to Characterize Parameter Sensitivity and Experiment to Study the Emphysema Progression in Clinical Images.

Parameter Sensitivity and Model Analysis

The results of the Exposure experiment are shown in Figures 5A,B. Figure 5A illustrates the effect of changing the number of particles inhaled for each simulated smoking exposure and the total time spent smoking. When exposure is zero, the model is able to capture that the tissue should remain healthy no matter how long the simulation runs. However, as both the exposure and total time spent smoking increase, the tissue starts getting damaged, independently on the values of the rate constants. For lower to medium exposures, the implicit stochasticity of the AB model and the variability of the rate constants lead to some the fluctuations of the results in function of the smoking time, but tissue life is always reduced by at least 50%. For higher exposure, tissue damage is irreversible and continues even after smoking cessation is simulated, as shown in Figure 5B. These outcomes nicely reflect the fact that smoke frequency and exposure are considered as one of the main risk factors for the development of COPD and emphysema (Yoshida and Tuder, 2007; Liu et al., 2008).

FIGURE 5

Figure 5. (A) Effect of changing the number of particles inhaled for each simulated smoking (exposure) and the total time of the simulation spent smoking (Smoking time). The value of each cell is the residual tissue life in % after the simulation stopped. (B) Emphysema progression in time. The parameters of the experiments are in figure (A): starting from experiment 0 in the lowest left corner to experiment 24 in the upper right.

In the model, the mechanical damage largely depends on the regulation of the collagen content, because the stiffness of this macromolecule is two orders of magnitude the stiffness of elastin. The degradation of collagen is heavily affected by TNFα through the recruitment of monocytes (Equation 1) and the activation of macrophages into M1 type (Equation 3) with positive feedback loops generated through IL1 (Equations 3, 7) and TNFα (Equation 6). In contrast, the activation of macrophages into M2 type is promoted by the anti-inflammatory cytokine IL10 (Equation 4), the production rate of which is positively retro-alimented by M2 macrophages (Equation 5). While IL10 inhibits TNF-alpha and IL1 (Equations 6, 7), it is also self-inhibited (Equation 5). Hence, the anti-inflammatory effect of IL10 is limited compared to the strong inflammatory effects of TNFα and IL1, because less positive feedback in favor of the promotion of type M2 activated macrophages. In the Parameters experiment, the exposure and total smoking time parameters were set to respectively 10 particles and 50 time steps to ensure that the system would be in a medium damage situation. Results were all very similar and the tissue life in each time step only varied with an average standard deviation of 0.166 units, revealing that the above interpretation of the model holds true regardless the variation of the rate constants within the considered ranges of values.

According to the analysis of the model equations, the promotion of the anabolic TGFβ (Equation 8) should be limited compared to the promotion of the catabolic MMP (Equation 12), and the persistent inflammation induced by particles should promote the unequivocate destruction of collagen (Equation 11). Nevertheless, we sometimes saw an increase in the mean amount of collagen. This outcome can be due to the mechanical feedback and mechanical tissue damage that promotes the secretion of TGF-beta and provides additional weight to collagen anabolism (Equations 9, 10). In our model emphysema progression, was indeed related to sustained inflammation that continued after smoking cessation (Willemse et al., 2005), but required the additional effect of DAMPS to relate inflammation, altered tissue turnover and tissue mechanics to cell endothelial death. This phenomenon needs to be further explored in a more mechanistic way, but our approximation of DAMPS effects allows qualitative validation of the simulated mechanisms for emphysema progression against clinical data (see below).

Emphysema Progression in Clinical Images

To test whether our model is able to produce images similar to those seen by clinicians, we use the public emphysema database described in section Experiment to Study the Emphysema Progression in Clinical Images. Images of some of the patches representative of the data we used in the experiment are presented in Figure 6. Parenchyma destruction in emphysema is strongly associated with decreased HU absorption value in CT images, and many image descriptors are commonly used to (semi-)automatically detect emphysema progression in CT images (Stern and Frank, 1994; Gevenois and Yernault, 1995; Madani et al., 2006). In the present study, emphysema progression is quantified through the well-known Mean Lung Density (MLD) (Heremans et al., 1992).

FIGURE 6

Figure 6. Database patches samples for the three categories of low (up), medium (middle) and severe (down) emphysema affectation. We see how the affectation is related to the appearance of bigger cluster of low attenuation areas from top to bottom.

We quantify all the patches from the database and group by the different degrees of emphysema severity. As it can be seen from Figure 7, images with increasing emphysema severity have also a lower MLD score. Differences of more than 40–60 HU between groups (1) and (2), (3) are significant with p-value of less than 0.01.

FIGURE 7

Figure 7. Relation between Mean Lung Density and emphysema progression. We used Mean Lung Density (MLD) to quantify patches belonging to the three emphysema classification levels. The figure shows that emphysema progression is associated with a mean lowering of the MLD values, due to the destruction of parenchyma and the diminishing of the CT attenuation value.

Once the association between emphysema progression and MLD score is determined, we take all the 69 patches annotated with low or no emphysema affectation and use them as input for our model. Images are quantified with MLD before and after model execution and the results are tested with t-test for statistical significance of the differences.

As we can see in Figure 8, there is a statistically significant difference of about 30 HU between the baseline and progression groups with a p-value of less than 0.001. We thus conclude that our implemented model is able to simulate changes that are in agreement with the progression of emphysema in clinical images quantified by MLD.

FIGURE 8

Figure 8. Effect of model execution. We run our emphysema progression model on 69 patches with low emphysema until completion. As expected, the resulting patches have a lower MLD than the original one, showing parenchyma destruction during model run. ***p-value < 0.001.

Scaling

The parallelization strategy for the AB part consists in assigning a job for each patch, as they were completely independent from each other, then recursively create a new job per voxel, which is the smaller unit we can parallelize for now. With 69 patches to process and 200 seeds per patch plus their four closest neighbors, this resulted in 69,000 jobs. The AB code uses OpenMP to parallelize the execution of the agents. For the FE part, we use the MPI capability of Elmer solver to partition the mesh and distribute the computation on 5 MPI processes per job. Finally, the coupling between AB and FE is executed in a sequential way with a python script. Most of the code executions for the coupling use libraries for which C code bindings were available (numpy, scipy, and skimage) and, thus, simulations run at almost native speed. During a single job execution, the computational time is taken mostly by the AB model (54%), then by the FE solver (37%) and finally by the coupling part (9%). Each job took about 40 min on the cluster and 1 h on a workstation computer (Intel i7 with 32 GB of RAM). A significant amount of time of the coupled simulations was spent on writing the files to disk to share data between the solvers. This could be reduced in future works by using faster SSD disk or in-memory access. All jobs would have taken months to be processed sequentially on the abovementioned workstation, but required 8 days on our cluster, by using a maximum of 256 simultaneous jobs. The use of the HPC resulted, therefore, in several orders of magnitude of computational time reduction, making the present study actually feasible. All jobs in the cluster use Sun Grid engine.

Conclusions and Future Works

In this paper we conceived, developed and tested a high performance multi-scale agent-based model of lung parenchyma evolution after repeated exposure to solid irritants such as the particles that arrive to the lung while smoking. We modeled the simplified behavior of immune system cells such as alveolar macrophages and neutrophils, and also cells in charge of wound healing mechanisms such as fibroblasts. Finally, the tissue behavior under the forces present in the lung during respiration was modeled using a FE elastic model. An initial analysis of sensitivity of the model to parameter variations confirmed (i) the ability of the model to point out particle inhalations as a major risk factor in emphysema pathogenesis, and (ii) the strong inertia of the catabolic shift of cell activity due to sustained inflammation that resulted in sustained damage to most of the tissue. A preliminary validation of the capacity of the model to cause a significant change was performed against clinical images on 69 cases of a public database of CT images affected by emphysema progression.

To the best of our knowledge, this model advances the state of the art because: (1) it includes a more detailed molecular model of inflammation and tissue remodeling (2) uses a FE solver to calculate the response to mechanical solicitations thus allowing for future extensions where arbitrary complex tissue constitutive equations could be used. (3) has a bi-directional coupling between AB and FE models (4) exploits HPC technologies so that enough tissue can be simulated to start validating against imaging data (5) uses clinical CT images to perform an initial validation of the capacity of the model. The implementation of a system of coupled ODEs into AB has the great advantage over a well-mixed model to take into account the spatial aspect, and the formation of self-sustaining spatial patterns that affect substantially the equilibrium points of the system (Brown et al., 2011).

The present model has, of course, several limitations. Simplifications were still made in the immune response and the mechanical model. In particular, the relative importance of DAMPS in the validated model suggest that more mechanistic development of this biological phenomenon are necessary. Additionally, while in this implementation of the model we used a 2D mapping between the alveolar exchange surface and the computational grid, in following works we will explore the effect of extending the connectivity of the tissue to 3D. On top of that, we do not account for heterogeneous tissue structures such as airways or blood vessels. However, we plan to do so in a following extension as the relevant information is already present in the CT images used to initialize the model. Effect of the mesh size and topology should be further explored. In a follow-up study we plan to automatically find the best parameters and better characterize the impact of the mesh on the stability of the solution with a convergence study.

Also, the validation is still somehow limited, as no histological comparison with ex-vivo animal models could be performed and the one on CT clinical images is limited to one clinical descriptor, namely the MLD score. As a future work, we are planning a retrospective study with COPD patients with 1-year follow-up. Of course, the real ground truth should be histology, which is unfortunately very difficult to obtain in human subjects. A promising alternative is to use mice models of emphysema.

We suggest that such a model, once properly extended and calibrated with histological and clinical data, could be useful to improve patient classification and prediction of exacerbations and thus contribute to the selection of a personalized therapy.

Author Contributions

MC: conceived the research, designed the model, generated and analyzed the experimental data, validated the implementation and wrote the paper; AO: implemented part of the ABM rules and helped with the collection and analysis of the experimental data; JN: helped with the preparation of the FE model with the writing of the paper and the analysis of the results; MG: oversaw the project and revised the paper. All authors read and approved the final version of the manuscript.

Funding

We gratefully acknowledge the NVIDIA Corporation for the donation of part of the hardware used for this research. This work was supported by the Spanish Ministry of Economy and Competitiveness (Project INSPIRE FIS2017-89535-C2-2-R, Ramon y Cajal contract RYC-2015-18888, Maria de Maeztu Units of Excellence Program MDM-2015-0502), and by the QUAES Foundation (Chair QUAES-UPF Computational Technologies for Healthcare).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

1. ^INSPIRE - Personalized computational models of COPD progression for patient phenotyping. FIS2017-89535-C2-2-R

2. ^Open source: gmsh.info

3. ^https://bitbucket.org/mrceresa/pandora

References

Abar, S., Theodoropoulos, G. K., Lemarinier, P., and O'Hare, G. M. P. (2017). Agent based modelling and simulation tools: a review of the state-of-art software. Comput. Sci. Rev. 24, 13–33. doi: 10.1016/j.cosrev.2017.03.001

CrossRef Full Text | Google Scholar

Akinbami, L. J., Moorman, J. E., Bailey, C., Zahran, H. S., King, M., Johnson, C. A., et al. (2012). Trends in asthma prevalence, health care use, and mortality in the United States 2001–2010. NCHS Data Brief 94, 1–8. Available online at: https://www.aap.org/en-us/Documents/medicalhome_resources_asthma_prevalence.pdf

Google Scholar

An, G. (2008). Introduction of an agent-based multi-scale modular architecture for dynamic knowledge representation of acute inflammation. Theor. Biol. Med. Modell. 5:11. doi: 10.1186/1742-4682-5-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Auffray, C., Fogg, D., Garfa, M., Elain, G., Join-Lambert, O., Kayal, S., et al. (2007). Monitoring of blood vessels and tissues by a population of monocytes with patrolling behavior. Science 317, 666–670. doi: 10.1126/science.1142883

PubMed Abstract | CrossRef Full Text | Google Scholar

Bellingan, G. J., Caldwell, H., Howie, S. E., Dransfield, I., and Haslett, C. (1996). In vivo fate of the inflammatory macrophage during the resolution of inflammation: inflammatory macrophages do not die locally, but emigrate to the draining lymph nodes. J. Immunol. 157, 2577–2585.

Google Scholar

Berger, L., Bordas, R., Burrowes, K., Grau, V., Tavener, S., and Kay, D. (2016). A poroelastic model coupled to a fluid network with applications in lung modelling. Int. J. Num. Methods. Biomed. Eng. 32:e02731. doi: 10.1002/cnm.2731

PubMed Abstract | CrossRef Full Text | Google Scholar

Birzle, A. M., Martin, C., Yoshihara, L., Uhlig, S., and Wall, W. A. (2018). Experimental characterization and model identification of the nonlinear compressible material behavior of lung parenchyma. J. Mechan. Behav. Biomed. Mater. 77, 754–763. doi: 10.1016/j.jmbbm.2017.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Bordas, R., Lefevre, C., Veeckmans, B., Pitt-Francis, J., Fetita, C., Brightling, C. E., et al. (2015). Development and analysis of patient-based complete conducting airways models. PLoS ONE 10:e0144105. doi: 10.1371/journal.pone.0144105

PubMed Abstract | CrossRef Full Text | Google Scholar

Brightling, C. E., Gupta, S., Gonem, S., and Siddiqui, S. (2012). Lung damage and airway remodelling in severe asthma. Clin. Exp. Allergy 42, 638–649. doi: 10.1111/j.1365-2222.2011.03917.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, B., Price, I. M., Toapanta, F. R., DeAlmeida, D. R., Wiley, C. A., Ross, T. M., et al. (2011) An agent-based model of inflammation fibrosis following particulate exposure in the lung. Math. Biosci. 231, 186–196.

Google Scholar

Burrowes, K. S., De Backer, J., Smallwood, R., Sterk, P. J., Gut, I., Wirix-Speetjens, R., et al. (2013). Multi-scale computational models of the airways to unravel the pathophysiological mechanisms in asthma and chronic obstructive pulmonary disease (AirPROM). Interface Focus 3:20120057. doi: 10.1098/rsfs.2012.0057

PubMed Abstract | CrossRef Full Text | Google Scholar

Calmet, H., Gambaruto, A. M., Bates, A. J., Vázquez, M., Houzeaux, G., and Doorly, D. J. (2016). Large-scale CFD simulations of the transitional and turbulent regime for the large human airways during rapid inhalation. Comput. Biol. Med.69, 166–180. doi: 10.1016/j.compbiomed.2015.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Castro, M., Rubin, A. S., Laviolette, M., Fiterman, J., De Andrade Lima, M., Shah, P. L., et al. (2010). Effectiveness and safety of bronchial thermoplasty in the treatment of severe asthma: a multicenter, randomized, double-blind, sham-controlled clinical trial. Am. J. Respir. Crit. Care Med. 181, 116–124. doi: 10.1164/rccm.200903-0354OC.

PubMed Abstract | CrossRef Full Text | Google Scholar

Ceresa, M., Olivares, A. L., Suelves, S. F., and González Ballester, M. A. (2017). “Multi-scale immunological and biomechanical model of emphysema progression,” in IEEE Proceedings of the 39th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (Seogwipo).

Google Scholar

Chernyavsky, I. L., Croisier, H., Chapman, L. A. C., Kimpton, L. S., Hiorns, J. E., Brook, B. S., et al. (2014). The role of inflammation resolution speed in airway smooth muscle mass accumulation in asthma: insight from a theoretical model. PLoS ONE 9:e90162. doi: 10.1371/journal.pone.0090162

PubMed Abstract | CrossRef Full Text | Google Scholar

Clancy, C. E., An, G., Cannon, W. R., Liu, Y., May, E. E., Ortoleva, P., et al. (2016). Multiscale modeling in the clinic: drug design and development. Ann. Biomed. Eng. 44, 2591–2610. doi: 10.1007/s10439-016-1563-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Cobbold, C. A., and Sherratt, J. A. (2000). Mathematical modelling of nitric oxide activity in wound healing can explain keloid and hypertrophic scarring. J. Theor. Biol. 204, 257–288. doi: 10.1006/jtbi.2000.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Cooper, B. G., Stockley, R. A., and Sapey, E. (2017). Small airways disease : time for a revisit? Int. J. Chron. Obstruct. Pulmon. Dis. 12, 2343–2353. doi: 10.2147/COPD.S138540

PubMed Abstract | CrossRef Full Text | Google Scholar

Cosio, M. G., Saetta, M., and Agusti, A. (2009). Immunologic aspects of chronic obstructive pulmonary disease. N. Eng. J. Med. 360, 2445–2454. doi: 10.1056/NEJMra0804752

PubMed Abstract | CrossRef Full Text | Google Scholar

Cristini, V., Koay, E. J., and Wang, Z. (2017). An Introduction to Physical Oncology: How Mechanistic Mathematical Modeling Can Improve Cancer Therapy Outcomes. CRC Press, Taylor & Francis Group. Available online at: https://www.crcpress.com/An-Introduction-to-Physical-Oncology-How-Mechanistic-Mathematical-Modeling/Cristini-Koay-Wang/p/book/9781466551343

Google Scholar

Crosby, L. M., and Waters, C. M. (2010). Epithelial repair mechanisms in the lung. AJP Lung Cell. Mol. Physiol. 298, L715–L731. doi: 10.1152/ajplung.00361.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Darby, I. A., Bisucci, T., Hewitson, T. D., and MacLellan, D. G. (1997). Apoptosis is increased in a model of diabetes-impaired wound healing in genetically diabetic mice. Int. J. Biochem. Cell Biol. 29, 191–200. doi: 10.1016/S1357-2725(96)00131-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Decramer, M., Janssens, W., and Miravitlles, M. (2012). Chronic obstructive lung disease. Lancet 379, 1341–1351. doi: 10.1016/S0140-6736(11)60968-9

CrossRef Full Text

Deisboeck, T. S., Zhang, L., Yoon, J., and Costa, J. (2009). In silico cancer modeling: is it ready for prime time?. Nat. Rev. Clin. Oncol. 6, 34–42. doi: 10.1038/ncponc1237

PubMed Abstract | CrossRef Full Text | Google Scholar

Deisboeck, T. S., Wang, Z., Macklin, P., and Cristini, V. (2011). Multiscale cancer modeling. Annu. Rev. Biomed. Eng. 13, 127–155. doi: 10.1146/annurev-bioeng-071910-124729

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Ventura, B., Lemerle, C., Michalodimitrakis, K., and Serrano, L. (2006). From in vivo to in silico biology and back. Nature 443, 527–533. doi: 10.1038/nature05127

PubMed Abstract | CrossRef Full Text | Google Scholar

Depalle, B., Qin, Z., Shefelbine, S. J., and Buehler, M. J. (2015). Influence of cross-link structure, density and mechanical properties in the mesoscale deformation mechanisms of collagen fibrils. J. Mech. Behav. Biomed. Mater. 52, 1–13. doi: 10.1016/j.jmbbm.2014.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Eberhardt, W., Akool, E. L. S., Rebhan, J., Frank, S., Beck, K. F., Franzen, R., et al. (2002). Inhibition of cytokine-induced matrix metalloproteinase 9 expression by peroxisome proliferator-activated receptor alpha agonists is indirect and due to a NO-mediated reduction of mRNA Stability. J. Biol. Chem.. 277, 33518–33528. doi: 10.1074/jbc.M202008200

PubMed Abstract | CrossRef Full Text | Google Scholar

Edwards, J. P., Zhang, X., Frauwirth, K. A., and Mosser, D. M. (2006). Biochemical and functional characterization of three activated macrophage populations. J. Leukoc. Biol. 80, 1298–1307. doi: 10.1189/jlb.0406249

PubMed Abstract | CrossRef Full Text | Google Scholar

Faner, R., and Agustí, Á. (2016). Multilevel, dynamic chronic obstructive pulmonary disease heterogeneity: a challenge for personalized medicine. Ann. Am. Thorac. Soc. 13, S466–S470. doi: 10.1513/AnnalsATS.201605-372AW

PubMed Abstract | CrossRef Full Text | Google Scholar

Folcik, V. A., An, G. C., and Orosz, C. G. (2007). The basic immune simulator: an agent-based model to study the interactions between innate and adaptive immunity. Theor. Biol. Med. Model. 4:39. doi: 10.1186/1742-4682-4-39

PubMed Abstract | CrossRef Full Text | Google Scholar

Freitas, R. K., and Schröder, W. (2008). Numerical investigation of the three-dimensional flow in a human lung model. J. Biomech. 41, 2446–2457. doi: 10.1016/j.jbiomech.2008.05.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Gershon, A. S., Warner, L., Cascagnette, P., Victor, J. C., and To, T. (2011). Lifetime risk of developing chronic obstructive pulmonary disease: a longitudinal population study. Lancet 378, 991–996. doi: 10.1016/S0140-6736(11)60990-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Gevenois, P. A., and Yernault, J. C. (1995). Can computed tomography quantify pulmonary emphysema?. Eur. Respir. J. 8, 843–848.

PubMed Abstract | Google Scholar

Goldklang, M., and Stockley, R. (2016). Pathophysiology of emphysema and implications. Chronic Obstr. Pulm. Dis. 3, 454–458. doi: 10.15326/jcopdf.3.1.2015.0175

PubMed Abstract | CrossRef Full Text | Google Scholar

Hehenberger, K., Heilborn, J. D., Brismar, K., and Hansson, A. (1998). Inhibited proliferation of fibroblasts derived from chronic diabetic wounds and normal dermal fibroblasts treated with high glucose is associated with increased formation of L-lactate. Wound Repair Regen. 6, 135–141. doi: 10.1046/j.1524-475X.1998.60207.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Heremans, A., Verschakelen, J. A., Van fraeyenhoven, L., and Demedts, M. (1992). Measurement of lung density by means of quantitative CT scanning: a study of correlations with pulmonary function tests. Chest 102, 805–811. doi: 10.1378/chest.102.3.805

PubMed Abstract | CrossRef Full Text | Google Scholar

Holgate, S. (2012). Trials and tribulations in identifying new biologic treatments for asthma. Trends Immunol. 33, 238–246. doi: 10.1016/j.it.2012.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Horio, T., Nishikimi, T., Yoshihara, F., Nagaya, N., Matsuo, H., Takishita, S., et al. (1998). Production and secretion of adrenomedullin in cultured rat cardiac myocytes and nonmyocytes: stimulation by interleukin-1α and tumor necrosis factor-α Endocrinology 139, 4576–4580.

Google Scholar

Huang, M., Sharma, S., Zhu, L., Keane, M., Luo, J., Zhang, L., et al. (2002). IL-7 inhibits fibroblast TGF-beta production and signaling in pulmonary fibrosis. J. Clin. Invest. 109, 931–937. doi: 10.1172/JCI0214685

PubMed Abstract | CrossRef Full Text | Google Scholar

Ignotz, R. A., and Massagué, J. (1986). Transforming growth factor-beta stimulates the expression of fibronectin and collagen and their incorporation into the extracellular matrix. J. Biol. Chem. 261, 4337–4345.

PubMed Abstract | Google Scholar

Jin, Y. F., Han, H. C., Berger, J., Dai, Q., and Lindsey, M. L. (2011). Combining experimental and mathematical modeling to reveal mechanisms of macrophage-dependent left ventricular remodeling. BMC Syst. Biol. 5:60. doi: 10.1186/1752-0509-5-60

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Y., and Lindsey, M. L. (2010). “Multi-scale modeling and analysis of left ventricular remodeling post myocardial infarction: integration of experimental and computational approaches,” in Application of Machine Learning (InTech), 267–280.

Google Scholar

Larsson, K. (2007). Aspects on pathophysiological mechanisms in COPD. J. Intern. Med. 262, 311–340. doi: 10.1111/j.1365-2796.2007.01837.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lenga, Y., Koh, A., Perera, A. S., McCulloch, C. A., Sodek, J., and Zohar, R. (2008). Osteopontin expression is required for myofibroblast differentiation. Circ. Res. 102, 319–327. doi: 10.1161/CIRCRESAHA.107.160408

PubMed Abstract | CrossRef Full Text | Google Scholar

Lintermann, A., and Schröder, W. (2017). Simulation of aerosol particle deposition in the upper human tracheobronchial tract. Eur. J. Mech. B/Fluids 63, 73–89. doi: 10.1016/j.euromechflu.2017.01.008

CrossRef Full Text | Google Scholar

Liu, Y., Lee, K., Perez-Padilla, R., Hudson, N. L., and Mannino, D. M. (2008). Outdoor and indoor air pollution and COPD-related diseases in high-and low-income countries [State of the Art Series. Chronic obstructive pulmonary disease in high-and low-income countries. Edited by G. Marks and M. Chan-Yeung. Number 2 in the series]. Int. J. Tuberc. Lung Dis. 12, 115–127.

PubMed Abstract | Google Scholar

Lotze, M. T., Zeh, H. J., Rubartelli, A., Sparvero, L. J., Amoscato, A. A., Washburn, N. R., et al. (2007). The grateful dead: damage-associated molecular pattern molecules and reduction/oxidation regulate immunity. Immunol. Rev. 220, 60–81. doi: 10.1111/j.1600-065X.2007.00579.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Maass, D. L., Hybki, D. P., White, J., and Horton, J. W. (2002). The time course of cardiac NF-[kappa]B activation and TNF-[alpha] secretion by cardiac myocytes after burn injury: contribution to burn-related cardiac contractile dysfunction. Shock. 17, 293–299. doi: 10.1097/00024382-200204000-00009

PubMed Abstract | CrossRef Full Text | Google Scholar

Madani, A., Zanen, J., De Maertelaer, V., and Gevenois, P. A. (2006). Pulmonary emphysema: objective quantification at multi–detector row CT—comparison with macroscopic and microscopic morphometry. Radiology 238, 1036–1043. doi: 10.1148/radiol.2382042196

PubMed Abstract | CrossRef Full Text | Google Scholar

Mantovani, A., Sica, A., Sozzani, S., Allavena, P., Vecchi, A., and Locati, M. (2004). The chemokine system in diverse forms of macrophage activation and polarization. Trends Immunol. 25, 677–686. doi: 10.1016/j.it.2004.09.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Marino, S., Hogue, I. B., Ray, C. J., and Kirschner, D. E. (2008). A methodology for performing global uncertainty and sensitivity analysis in systems biology. J. Theor. Biol. 254, 178–196. doi: 10.1016/j.jtbi.2008.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinez, C. H., Chen, Y. H., Westgate, P. M., Liu, L. X., Murray, S., Curtis, J. L., et al. (2012). Relationship between quantitative CT metrics and health status and BODE in chronic obstructive pulmonary disease. Thorax 67, 399–406. doi: 10.1136/thoraxjnl-2011-201185

PubMed Abstract | CrossRef Full Text | Google Scholar

Matzinger, P. (2002). The danger model: a renewed sense of self. Science 296, 301–305. doi: 10.1126/science.1071059

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, F., and Lowell, C. A. (1997). Lipopolysaccharide (LPS)-induced macrophage activation and signal transduction in the absence of Src-Family Kinases Hck, Fgr, and Lyn. J. Exp. Med. 185, 1661–1670. doi: 10.1084/jem.185.9.1661

PubMed Abstract | CrossRef Full Text | Google Scholar

Mishima, M., Hirai, T., Itoh, H., Nakano, Y., Sakai, H., Muro, S., et al. (1999). Complexity of terminal airspace geometry assessed by lung computed tomography in normal subjects and patients with chronic obstructive pulmonary disease. Proc. Natl. Acad. Sci. U. S. A. 96, 8829–8834.

PubMed Abstract | Google Scholar

Mizuno, S., Ishizaki, T., Kadowaki, M., Akai, M., Shiozaki, K., Iguchi, M., et al. (2017). p53 signaling pathway polymorphisms associated with emphysematous changes in patients with COPD. Chest 152, 58–69. doi: 10.1016/j.chest.2017.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Moro, L., Arbini, A. A., Yao, J. L., di Sant'Agnese, P. A., Marra, E., and Greco, M. (2008). Loss of BRCA2 promotes prostate cancer cell invasion through up-regulation of matrix metalloproteinase-9. Cancer Sci. 99, 553–563. doi: 10.1111/j.1349-7006.2007.00719.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ochs, M., Nyengaard, J. R., Jung, A., Knudsen, L., Voigt, M., Wahlers, T., et al. (2004). The number of alveoli in the human lung. Am. J. Respir. Crit. Care Med. 169, 120–124. doi: 10.1164/rccm.200308-1107OC

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliver, J. C., Bland, L., Oettinger, C., Arduino, M., McAllister, S., Aguero, S., et al. (1993). Cytokine kinetics in an in vitro whole blood model following an endotoxin challenge. Lymphokine Cytokine Res. 12:115.

PubMed Abstract | Google Scholar

Onozaki, K., Urawa, H., Tamatani, T., Iwamura, Y., Hashimoto, T., Baba, T., et al. (1988). Synergistic interactions of interleukin 1, interferon-beta, and tumor necrosis factor in terminally differentiating a mouse myeloid leukemic cell line (M1). Evidence that interferon-beta is an autocrine differentiating factor. J. Immunol. 140, 112–119

PubMed Abstract | Google Scholar

Parameswaran, H., Majumdar, A., and Suki, B. (2011). Linking microscopic spatial patterns of tissue destruction in emphysema to macroscopic decline in stiffness using a 3D computational model. PLoS Comput. Biol. 7:e1001125. doi: 10.1371/journal.pcbi.1001125

PubMed Abstract | CrossRef Full Text | Google Scholar

Partridge, M. R., van der Molen, T., Myrseth, S. E., and Busse, W. W. (2006). Attitudes and actions of asthma patients on regular maintenance therapy: the INSPIRE study. BMC Pulm. Med. 6:13. doi: 10.1186/1471-2466-6-13

PubMed Abstract | CrossRef Full Text | Google Scholar

Pavord, I. D., Korn, S., Howarth, P., Bleecker, E. R., Buhl, R., Keene, O. N., et al. (2012). Mepolizumab for severe eosinophilic asthma (DREAM): a multicentre, double-blind, placebo-controlled trial. Lancet 380, 651–659. doi: 10.1016/S0140-6736(12)60988-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Porcheray, F., Viaud, S., Rimaniol, A. C., Léone, C., Samah, B., Dereuddre-Bosquet, N., et al. (2005). Macrophage activation switching: an asset for the resolution of inflammation. Clin. Exp. Immunol. 142, 481–489 doi: 10.1111/j.1365-2249.2005.02934.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Råback, P., Malinen, M., Ruokolainen, J., Pursula, A., and Zwinger, T. (2013). Elmer Models Manual. Helsinki: CSC-IT Center for Science.

Google Scholar

Repapi, E. (2010). Genome-wide association study identifies five loci associated with lung function. Nat. Genet. 42, 36–44. doi: 10.1038/ng.501

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, C. J., Yoshihara, L., Ismail, M., and Wall, W. A. (2017a). Computational modelling of the respiratory system: discussion of coupled modelling approaches and two recent extensions. Comput. Methods Appl. Mech. Eng. 314, 473–493. doi: 10.1016/j.cma.2016.08.010

CrossRef Full Text | Google Scholar

Roth, C. J., Yoshihara, L., and Wall, W. A. (2017b). A simplified parametrised model for lung microstructures capable of mimicking realistic geometrical and mechanical properties. Comput. Biol. Med. 89, 104–114. doi: 10.1016/j.compbiomed.2017.07.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubio-Campillo, X. (2014). “Pandora: a versatile agent-based modelling platform for social simulation,” in Proceedings of SIMUL 2014, The Sixth International Conference on Advances in System Simulation (Nice), 29–34.

Google Scholar

Sharafkhaneh, A., Hanania, N. A., and Kim, V. (2008). Pathogenesis of emphysema: from the bench to the bedside. Proc. Am. Thorac. Soc. 5, 475–477. doi: 10.1513/pats.200708-126ET

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherman, V. R., Yang, W., and Meyers, M. A. (2015). The materials science of collagen. J. Mech. Behav. Biomed. Mater. 52, 22–50. doi: 10.1016/j.jmbbm.2015.05.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Sorensen, L., Shaker, S. B., and De Bruijne, M. (2010). Quantitative analysis of pulmonary emphysema using local binary patterns. IEEE Trans. Med. Imaging 29, 559–569. doi: 10.1109/TMI.2009.2038575

PubMed Abstract | CrossRef Full Text | Google Scholar

Stefan, M. I., and Le Novère, N. (2013). Cooperative binding. PLoS Comput. Biol. 9:e1003106. doi: 10.1371/journal.pcbi.1003106

PubMed Abstract | CrossRef Full Text | Google Scholar

Steinmüller, C., Franke-Ullmann, G., Lohmann-Matthes, M. L., and Emmendorffer, A. (2000). Local activation of nonspecific defense against a respiratory model infection by application of interferon-gamma. Comparison between rat alveolar and interstitial lung macrophages. Am. J. Respir. Cell Mol. Biol. 22, 481–490. doi: 10.1165/ajrcmb.22.4.3336

PubMed Abstract | CrossRef Full Text | Google Scholar

Stern, E. J., and Frank, M. S. (1994). CT of the lung in patients with pulmonary emphysema: diagnosis, quantification, and correlation with pathologic and physiologic findings. AJR. Am. J. Roentgenol. 162, 791–798. doi: 10.2214/ajr.162.4.8140992

PubMed Abstract | CrossRef Full Text | Google Scholar

Suki, B., Lutchen, K. R., and Ingenito, E. P. (2003). On the progressive nature of emphysema: roles of proteases, inflammation, and mechanical forces. Am. J. Resp. Crit. Care Med. 168, 516–521. doi: 10.1164/rccm.200208-908PP

PubMed Abstract | CrossRef Full Text | Google Scholar

Suki, B., Stamenović, D., and Hubmayr, R. (2011). Lung parenchymal mechanics. Compr. Physiol. 1, 1317–1351 doi: 10.1002/cphy.c100033

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanaka, N., Hoshino, Y., Gold, J., Hoshino, S., Martiniuk, F., Kurata, T., et al. (2005). Interleukin-10 induces inhibitory C/EBβ through STAT-3 and represses HIV-1 transcription in macrophages. Am. J. Respir. Cell Mol. Biol. 33, 406–411. doi: 10.1165/rcmb.2005-0140OC

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsutsumi, Y., Kihira, T., Tsunoda, S., Kamada, H., Nakagawa, S., Kaneda, Y., et al. (1996). Molecular design of hybrid tumor necrosis factor-alpha III: polyethylene glycol-modified tumor necrosis factor-alpha has markedly enhanced antitumor potency due to longer plasma half-life and higher tumor accumulation. J. Pharmacol. Exp. Ther. 278:1006

PubMed Abstract | Google Scholar

Verdugo, F., Roth, C. J., Yoshihara, L., and Wall, W. A. (2017). Efficient solvers for coupled models in respiratory mechanics. Int. J. Numer. Method. Biomed. Eng. 33:e02795. doi: 10.1002/cnm.2795

PubMed Abstract | CrossRef Full Text | Google Scholar

Vestbo, J., Hurd, S. S., Agustí, A. G., Jones, P. W., Vogelmeier, C., Anzueto, A., et al. (2013). Global strategy for the diagnosis, management, and prevention of chronic obstructive pulmonary disease GOLD executive summary. Am. J. Respir. Crit. Care Med. 187, 347–365. doi: 10.1164/rccm.201204-0596PP

PubMed Abstract | CrossRef Full Text | Google Scholar

Wahl, S. M., Hunt, D. A., Wakefield, L. M., McCartney-Francis, N., Wahl, L. M., Roberts, A. B., et al. (1987). Transforming growth factor type beta induces monocyte chemotaxis and growth factor production. Proc. Natl. Acad. Sci.U.S.A. 84, 5788–5792. doi: 10.1073/pnas.84.16.5788

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Yang, T., Ma, Y., Halade, G. V., Zhang, J., Lindsey, M. L., et al. (2012). Mathematical modeling and stability analysis of macrophage activation in left ventricular remodeling post-myocardial infarction. BMC Genomics 13(Suppl. 6):S21. doi: 10.1186/1471-2164-13-S6-S21

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., and Maini, P. K. (2017). Editorial: special section on multiscale cancer modeling. IEEE Trans. Biomed. Eng. 64, 501–503. doi: 10.1109/TBME.2017.2655439

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Butner, J. D., Kerketta, R., Cristini, V., and Deisboeck, T. S. (2015). Simulating cancer growth with multiscale agent-based modeling. Semin. Cancer Biol. 30, 70–78. doi: 10.1016/j.semcancer.2014.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Wiechert, L., Comerford, A., Rausch, S., and Wall, W. A. (2011). “Advanced multi-scale modelling of the respiratory system,” in Fundamental Medical and Engineering Investigations on Protective Artificial Respiration (Berlin; Heidelberg: Springer-Verlag), 1–32.

Google Scholar

Willemse, B. W., ten Hacken, N. H., Rutgers, B., Lesman-Leegte, I. G., Postma, D. S., and Timens, W. (2005). Effect of 1-year smoking cessation on airway inflammation in COPD and asymptomatic smokers. Eur. Respir. J. 26, 835–845. doi: 10.1183/09031936.05.00108904

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshida, T., and Tuder, R. M. (2007). Pathobiology of cigarette smoke-induced chronic obstructive pulmonary disease. Physiol. Rev. 87, 1047–1082. doi: 10.1152/physrev.00048.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Yousem, S. A. (2006). Respiratory bronchiolitis-associated interstitial lung disease with fibrosis is a lesion distinct from fibrotic nonspecific interstitial pneumonia: A proposal. Mod. Pathol. 19, 1474–1479. doi: 10.1038/modpathol.3800671

PubMed Abstract | CrossRef Full Text | Google Scholar

Zemans, R. L., Henson, P. M., Henson, J. E., and Janssen, W. J. (2015). Conceptual approaches to lung injury and repair. Ann. Am. Thor. Soc. 12(Suppl. 1), S9–S15. doi: 10.1513/AnnalsATS.201408-402MG

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Ahmad, M., and Gronowicz, G. (2003). Effects of transforming growth factor-beta 1 (TGF-[beta]1) on in vitro mineralization of human osteoblasts on implant materials. Biomaterials 24, 2013–2020. doi: 10.1016/S0142-9612(02)00616-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, D., Zhou, Y., Jiang, C., Zhao, Z., He, F., and Ran, P. (2017). Small airway disease: a different phenotype of early stage COPD associated with biomass smoke exposure. Respirology 23, 198–205. doi: 10.1111/resp.13176

PubMed Abstract | CrossRef Full Text | Google Scholar