A computational model of invasive aspergillosis in the lung and the role of iron

Background Invasive aspergillosis is a severe infection of immunocompromised hosts, caused by the inhalation of the spores of the ubiquitous environmental molds of the Aspergillus genus. The innate immune response in this infection entails a series of complex and inter-related interactions between multiple recruited and resident cell populations with each other and with the fungal cell; in particular, iron is critical for fungal growth. Results A computational model of invasive aspergillosis is presented here; the model can be used as a rational hypothesis-generating tool to investigate host responses to this infection. Using a combination of laboratory data and published literature, an in silico model of a section of lung tissue was generated that includes an alveolar duct, adjacent capillaries, and surrounding lung parenchyma. The three-dimensional agent-based model integrates temporal events in fungal cells, epithelial cells, monocytes, and neutrophils after inhalation of spores with cellular dynamics at the tissue level, comprising part of the innate immune response. Iron levels in the blood and tissue play a key role in the fungus’ ability to grow, and the model includes iron recruitment and consumption by the different types of cells included. Parameter sensitivity analysis suggests the model is robust with respect to unvalidated parameters, and thus is a viable tool for an in silico investigation of invasive aspergillosis. Conclusions Using laboratory data from a mouse model of invasive aspergillosis in the context of transient neutropenia as validation, the model predicted qualitatively similar time course changes in fungal burden, monocyte and neutrophil populations, and tissue iron levels. This model lays the groundwork for a multi-scale dynamic mathematical model of the immune response to Aspergillus species. Electronic supplementary material The online version of this article (doi:10.1186/s12918-016-0275-2) contains supplementary material, which is available to authorized users.

Grid cells. The airway is a tube of grid cells which branches about one-third of the way across the horizontal span of the model. In total, approximately 16% of the grid cells make up the airway. Airway grid cells have no properties other than the number of spores at that location. Non-specific interstitial tissue cells comprise approximately 80% of the grid cells. The state variables for these cells are iron, macrophage-specific cytokines, and neutrophil-specific cytokines. Four blood vessels run the length of the model, made up of blood cells occupying one grid cell each; they comprise approximately 4% of the grid cells. These cells serve as the recruitment sites for macrophages and neutrophils as well as the source of iron for the interstitial space. Iron is introduced in the blood cells and diffuses into surrounding tissue. Cytokine levels are first affected by the epithelial cells and diffuse from the epithelium throughout the interstitial and blood cells. These cytokines serve as aggregations of macrophageand neutrophil-specific cytokines, such as TNF-α, IL-8, IFN γ , and others.
Fungal spores. A. fumigatus spores drift through the airway, and if a spore meets the epithelial cell wall it lodges there with a fixed probability. This probability is small, as ciliary beating serves as the primary mechanism for elimination of the fungus [4]. Once lodged, a spore cannot be swept away. Spores that are swept away are placed on the nearest airway grid cell and continue to drift. Once a spore has traversed the entire airway, it is removed from the simulation.
Approximately 30% of spores are internalized by an epithelial cell [5]. Since approximately 3% of internalized spores survive and 34% of these germinate [6], internalized spores germinate with approximate probability 0.01. Internalized or not, all lodged spores enter a resting stage. Upon completion of the resting stage, all non-internalized spores become swollen. After remaining in the swollen stage for some time, spores germinate. Since the resting stage lasts roughly 2 hours and germination occurs between 6 -8 hours [7], the swelling process lasts approximately 5 hours. A germinated spore spawns a fungal hyphal cell, which is connected to the spore and grows into neighboring tissue in a random direction. Since behavior of hyphae are different from that of fungal spores, these agents are described separately.
Epithelial cells. Epithelial cells form the boundary between the airway and the interstitial tissue. Each time step, these cells first determine the number of fungal spores and hyphae within their detection radius. These counts indicate the amount of macrophage-specific and neutrophil-specific cytokines to produce, respectively. While the cells produce the cytokine levels based on fungal presence, these levels are in fact grid cell state variables, not attributes of the epithelial cells. After adjusting local cytokine levels accordingly, each epithelial cell damages any internalized spores. This continues until either the spore germinates and kills the epithelial cell, or the spore dies. It is observed in [8] that only 3% of conidia survive at least 36 hours -hence 30 hours is chosen as the approximate time it takes an epithelial cell to kill an internalized conidium.
Macrophages. Macrophages are recruited to the site of infection via the bloodstream. Each time step, if there is at least one blood cell whose macrophage-specific cytokine level is above the macrophage recruitment threshold, there is a possibility of a macrophage spawning at exactly one such location. All macrophages absorb some of the macrophage cytokines on their current location, simulating uptake by receptors. Macrophages have a detection radius for conidial spores and hyphae; if a macrophage has fewer than two internalized spores, it will internalize a spore within this radius. Internalization prevents spores from growing hyphae. The internalization process takes approximately two hours [9]; since engulfed conidia are assumed unable to escape, this is used as the time it takes for a macrophage to kill an internalized spore. Macrophages do not phagocytose hyphae, but they are known to produce IL-8 to aid in neutrophil recruitment [10,11] -this is simulated by having macrophages increase the neutrophil-specific cytokine level at their current location based on the amount of nearby hyphae.
Fungal hyphae. Hyphae grow out of germinated fungal spores. They are attached to the spores and thus do not move. Once a hyphal cell appears, it spends time in a growth stage, during which no new hyphae can grow out of it. During the growth phase, if the hyphae has not been internalized, it attempts to absorb iron from surrounding tissue. Once the cell obtains enough iron, a new hyphal cell (or two new cells, if branching occurs) grows from the parent cell, with the new cells inheriting iron from the parent. Once it has spawned a new cell, a parent cell is no longer eligible for spawning future cells, though it does continue to take up iron from the environment. Many hyphae state variable values are unitless, abstract measures (e.g., iron levels), while others are chosen empirically.
Neutrophils. Neutrophils are recruited from the bloodstream in much the same way as macrophages -if there is at least one blood cell with neutrophil-specific cytokines above a certain threshold, then there is the possibility of a neutrophil spawning at exactly one such location. To indicate absorption of cytokines by receptors, every neutrophil reduces the neutrophil-specific cytokine level at its current location. As long as there are grid cells nearby with cytokine levels above the recruitment threshold, neutrophils move towards the grid cell with highest cytokine level -this simulates the chemotactic movement of recruited neutrophils [12]. Otherwise, they move to the center of a randomly chosen neighboring non-airway grid cell. Neutrophils raise cytokine levels in order to boost recruitment of additional neutrophils by an amount equal to the number of hyphae within their detection radius. Neutrophils are initialized with a fixed number of granules, which are deposited every time step in order to degrade and kill hyphae. Various studies indicate that the approximate killing time of hyphae by neutrophils is 2 hours [13,14]. All hyphae within the neutrophil detection radius are damaged as long as the neutrophil has not run out of granules; during this process iron is removed from the neutrophil's location in order to simulate sequestration. Additionally, neutrophils do not move while they are degrading nearby hyphae. Recruited neutrophils have a lifespan of 1 -2 days [15]; thus in the model, the average lifespan is set at 36 hours.

A.3 Process overview and scheduling
In the simulation, time is discrete. Entity routines are executed serially -that is, one entity executes every command in the routine, then the next entity does the same, and so on. This means that entity state variables are updated asynchronously; pertinent state variables are updated as each command is executed. The order in which entities perform the routine is randomly chosen at each time step. The pseudocode in Algorithm S3.1 describes the global process overview; sub-processes are provided in Section A.7. State variables referred to in the algorithms are described in Table S2.
Algorithm S3.1 Pseudocode for the entire model process.
1: setup map and constants 2: while current sim time < total sim time do 3: conidia routine: if mobile then 6: if at edge of model then die else conidia move end if

A.4 Design concepts
Basic principles. The biological mechanics of cell-cell interactions are the basic principles underlying the design of this model. These are described in literature and derived from experimentation, both in vivo and in vitro. Some behavioral processes, such as phagocytosis, occur on a local scale between two cells, while others, such as recruitment of macrophages and neutrophils, are tissue-level processes. The model aims to tie together what is known about individual cell behavior to create an in silico model of a larger-scale process.
Emergence. The ability of the fungus to survive under various conditions is an emergent property of individual fungal cells competing for iron. At the same time, the tendency for survival of infection is an emergent property of the entire immune response.
Adaptation. Macrophages and neutrophils both generally travel in the direction of highest cytokine concentration; thus, they adapt their movement in response to developing infection.
Objectives. The objective of immune cells is to remove all fungal cells, which is measured simply as a count of those cells. Thus whether a patient lives or dies depends on the amount of fungus present in the system.

Learning. Agents do not change adaptive traits over time in this simulation.
Prediction. Agent prediction is similar to adaptation in this case: macrophages and neutrophils inherently predict that following the path of highest cytokine concentration will lead to the area where immune response is most needed.
Sensing. Fungal spores sense the type of grid cell they are near, be it airway or interstitial space. In addition, both spores and hyphae sense the amount of iron present at their current location, and use this information to 'decide' if they will grow. Epithelial cells, macrophages, and neutrophils can all sense nearby spores and hyphae as well as the cytokine levels of all neighboring cells. All sensing is local.
Interaction. While section A.2 provides a fairly comprehensive overview of entity interactions, it is perhaps helpful to mention several interactions which do not occur. Specifically, macrophages and neutrophils do not directly interact with other immune cells, either of their type or another. For example, the presence of many macrophages at one location has no bearing on the likelihood of another macrophage targeting that location. Fungal hyphae interact with other hyphae as long as they are connected -in particular, iron uptake is shared among all hyphae at a given location, and when new hyphae are spawned, they inherit iron from the parent cell.
Stochasticity. Agent movement, when not dictated by concerns such as cytokine levels, is stochastic: agents face in the direction of a randomly chosen legal space (i.e., non-airway) and move in that direction. Hyphae grow in random directions, leading to clumps of fungus rather than long strands; hyphal branching is stochastic as well. In general, stochasticity is employed as a substitute for unknown mechanisms of the biological entities being represented.
Collectives. Fungal hyphae form a collective of sorts, as they maintain information about every other hyphae they are connected to. Iron stores are shared among local neighbors, which has an indirect effect on the collective group. These groups arise naturally and from the local rules of individual entities; they are not entities themselves and have no distinct state variables.
Observation. Many data can be collected from the model for analysis and understanding. These include fungal cell counts, iron levels, cytokine levels, macrophage and neutrophil counts, and collateral tissue damage. Any or all of these may be collected during each simulation, depending on the needs and interests of the researcher. The data are collated in a universal spreadsheet format for ease of statistical interpretation and analysis.

A.5 Initialization
Many initialization values are given in Table S2; as such, they are not repeated here. In addition to these values, the following information may be necessary to reproduce results: blood and interstitial tissue cells begin with iron at 20; all grid cells begin with cytokine levels at 0. The initial inoculum is placed at one end of the airway. There are no macrophages or neutrophils until recruitment occurs via cytokine triggering. Macrophages and neutrophils are created at the blood cell where recruitment occurs. Given that the model is meant to function as an in silico laboratory, it is not meant to be entirely robust with respect to initial conditions. Indeed, a key feature of the model is the ability to investigate the effect of initial conditions on outcome.

A.6 Input data
The model requires a map file (in .txt format) in order to set up grid cell types and locations. The map file is a three-dimensional matrix. Grid cell layout is determined aesthetically to serve the visual representation of model dynamics. No other input data is used.

A.7 Submodels
All submodels are presented here as pseudocode. The notation rand(a, b) indicates a uniformly distributed random number chosen from interval (a, b). for all non-captive hyphal or hyphal-tip fungus here with iron f < iron max (f ) do 3: iron f = iron f + (iron abs(f ) · iron p /no. fungus here) iron is split among all eligible fungus 4: iron f = min{iron f , iron max (f )} iron capped at iron max (f ) 5: end for 6: iron p = (1 − iron abs(f )) · iron p grid cell iron reduced 7: end for Algorithm S3.4 grow hyphae 1: for all germinated non-captive spores which haven't spawned hyphae do 2: spawn hyphal cell 3: hyphal cell routine: 4: move forward spacing in a random non-airway direction 5: set stage to 'hyphal tip' 6: set growth time to t grow begin swelling phase 7: iron f = iron f of parent cell 8: end routine 9: end for 10: for all hyphal tip cells with iron f > iron min (f ) and growth time reached do 11: num to spawn = 1 12: if rand(0, 1) < p branch then probability of branching 13: num to spawn = 2 14: iron p = min{(iron p + iph · time step/60), iron max (p)} iron is capped at iron max (p) 3: end for 4: for all non-airway grid cells do 5: give each neighbor cell (iron dif · iron p /26) iron 26 neighbors in 3D cyto m = cyto m + cyto rate · spore count increase cytokine level based on counts 5: cyto n = cyto n + cyto rate · (hyphae count + spore count) 6: end for Algorithm S3.7 epithelial damage internal spores 1: for all internalized spores do 2: add cyto m /26 to each neighboring cell cyto m refers to cytokine level of contributing cell 5: add cyto n /26 to each neighboring cell diffusion to 26 neighbors in 3D 6: end for 7: airway cells: cyto m , cyto n = 0 no cytokines in airway Algorithm S3.9 check for new macrophages 1: Repeat twice: 2: if any blood cells with cyto m ≥ m recr then 3: let L be one such location, selected at random 4: if rand(0, 1) < p recr (m) then probability of recruitment 5: spawn new macrophage at location L set num reps = (current sim time − 48)/8 suppress neutrophil production between days 2 and 4 4: end if 5: Repeat num reps: 6: if any blood cells with cyto n ≥ n recr then 7: let L be one such location, selected at random 8: spawn new neutrophil at location L 9: end if Algorithm S3.15 neutrophil absorb cytokines 1: cyto n = (1 − n abs ) · cyto n Algorithm S3.16 neutrophil produce neutrophil cytokine 1: let h be number of non-internalized fungus in stage 'hyphal tip' or 'hyphal' within n det 2: grid cell here: cyto n = cyto n + N n · h Algorithm S3.17 neutrophil move 1: if degranulating nearby hyphae then 2: do not move 3: else if any neighboring grid cells with cyto n ≥ n recr then 4: move to neighbor with highest cyto n follow highest cytokine concentration 5: else 6: move to random neighboring grid cell may wrap around map 7: end if Algorithm S3.18 neutrophil damage hyphae 1: if any hyphae within n det and gran > 0 then 2: hyphae: health f = health f − (init health(f ) · time step/n kill ) all grid cells within n det : set iron p = 0 neutrophils sequester iron 5: end if