ABSTRACT
Mycobacterium tuberculosis, the causative agent of tuberculosis (TB), induces formation of granulomas, structures in which immune cells and bacteria colocalize. Macrophages are among the most abundant cell types in granulomas and have been shown to serve as both critical bactericidal cells and targets for M. tuberculosis infection and proliferation throughout the course of infection. Very little is known about how these processes are regulated, what controls macrophage microenvironment-specific polarization and plasticity, or why some granulomas control bacteria and others permit bacterial dissemination. We take a computational-biology approach to investigate mechanisms that drive macrophage polarization, function, and bacterial control in granulomas. We define a “macrophage polarization ratio” as a metric to understand how cytokine signaling translates into polarization of single macrophages in a granuloma, which in turn modulates cellular functions, including antimicrobial activity and cytokine production. Ultimately, we extend this macrophage ratio to the tissue scale and define a “granuloma polarization ratio” describing mean polarization measures for entire granulomas. Here we coupled experimental data from nonhuman primate TB granulomas to our computational model, and we predict two novel and testable hypotheses regarding macrophage profiles in TB outcomes. First, the temporal dynamics of granuloma polarization ratios are predictive of granuloma outcome. Second, stable necrotic granulomas with low CFU counts and limited inflammation are characterized by short NF-κB signal activation intervals. These results suggest that the dynamics of NF-κB signaling is a viable therapeutic target to promote M1 polarization early during infection and to improve outcome.
INTRODUCTION
The pathological hallmark of the immune response to Mycobacterium tuberculosis is the formation of organized clusters of immune cells and mycobacteria called granulomas. Granulomas typically form in lungs of an infected host following inhalation of M. tuberculosis and the establishment of infection. These complex and dynamic structures can immunologically restrain mycobacterial proliferation and physically contain bacterial dissemination and thus represent an important determinant of disease outcome (1–3). Although there are several types of observed granulomas, the classical necrotic (caseous) granuloma has a necrotic core that may contain extracellular bacteria surrounded by an epithelioid macrophage-rich region, which may contain infected macrophages, and an outer cuff of lymphocytes intermixed with macrophages (4, 5). This unique spatial organization emerges as early as 1 month after infection in primates, with many granuloma lesions forming in the lung, each of a different size and bacterial burden (6, 7).
Macrophages are the most abundant cells in granulomas and have been shown to play a key role throughout the course of infection in all infection models, including humans and the nonhuman primate (8). Multiple other cell types, including neutrophils (9–13), dendritic cells (DCs), and fibroblasts, can also be found but are not the most abundant cell types (4). Macrophages perform many functions: they kill M. tuberculosis, serve as reservoirs of M. tuberculosis intracellular infection and proliferation, and regulate immune responses by producing a variety of pro- and anti-inflammatory cytokines (3, 8). Distinct macrophage activation states (polarization) have been classified as classical M1 activation (stimulated by Toll-like receptor [TLR] ligands and gamma interferon [IFN-γ], signaling via STAT1) or alternative M2 activation (stimulated by interleukin-4 [IL-4] and interleukin-13 [IL-13] via STAT6, as well as interleukin-10 [IL-10] via STAT3) (14–17). The M1, or proinflammatory, phenotype is characterized by expression of high levels of proinflammatory cytokines (e.g., tumor necrosis factor alpha [TNF-α]), high production of reactive nitrogen and oxygen intermediates, promotion of a Th1 response, and strong microbicidal and tumoricidal activity (15, 16, 18). Resistance to intracellular pathogens is generally considered to be due to the action of M1-like macrophages, which characterize the early phases of infection with M. tuberculosis in mice (19, 20). The M2 or anti-inflammatory phenotype is characterized by expression of anti-inflammatory cytokines (e.g., IL-10) and the promotion of tissue healing and remodeling, as well as by immunoregulatory functions (18, 20). Experimental studies have shown that these polarized macrophage phenotypes (M1 and M2) can be reversible both in vitro and in vivo (21–24). It has been reported that a switch between pro- and anti-inflammatory cytokine profiles can be observed during the transition from acute to chronic infection and may be a regulatory mechanism providing protection against excessive inflammation (25). On the other hand, pathogens that have evolved strategies to interfere with M1-associated killing can drive a phenotypic switch toward M2 polarization (26–29).
TNF and IFN-γ are key proinflammatory mediators that are elicited during M. tuberculosis infection. Both experimental and modeling studies, as well as clinical data from humans (i.e., anti-TNF therapies), have revealed that TNF plays a major role in host defense against M. tuberculosis, in both the active and chronic phases of infection, as well in proper granuloma function (30–35). A necessary step toward a protective immune response to mycobacterial infection requires NF-κB activation of macrophages (36, 37). However, many studies have also shown how specific NF-κB-mediated pathways may be exploited by bacterial pathogens to promote survival (36, 38, 39). NF-κB activation of macrophages occurs after cell surface ligation of TNF by tumor necrosis factor 1 (TNFR1) (40) and/or engagement of Toll-like receptors (TLRs) by microbial ligands (41–44). The signaling pathways downstream of NF-κB include the induction of a variety of proinflammatory genes along with upregulation of phagolysosome fusion-mediated killing of mycobacteria (45). Macrophage bactericidal activity and phagocytosis are significantly enhanced by IFN-γ when it acts in concert with TNF and have been shown to be necessary in both humans and mice for control of M. tuberculosis (32, 46–48). IFN-γ-inducible STAT1 phosphorylation transcriptionally upregulates production of inducible nitric oxide synthase (iNOS), which enables production of nitric oxide and related reactive nitrogen intermediates (RNI) by macrophages, a major antimycobacterial effector mechanism that is characteristic of proinflammatory M1-like macrophages (49, 50).
One of the primary downregulators of the immune response is IL-10 (4, 51–55); however, its role in antimycobacterial immunity remains to be fully elucidated. In humans, high levels of IL-10 are found in lungs, serum, sputum, and bronchoalveolar lavage (BAL) fluid of patients with pulmonary tuberculosis (TB) (4, 56). There are data suggesting a necessary role in control of TB (57–61), while others suggest that IL-10 plays no part in the dynamics (62–68). Our modeling studies predict a protective role for IL-10, balancing inflammation, limiting tissue damage, and regulating the dominant macrophage phenotype in the lymph node and lung (69–71). IL-10 induces STAT3 activation and downstream gene expression when it binds with the IL-10 receptor (IL-10R) complex, which is associated with M2-polarized macrophages. Furthermore, ligation of IL-10 with its receptor complex directly inhibits the TNF mRNA transcription pathway through STAT3-related factors (72, 73).
Overall, macrophage function, diversity, and location within a granuloma environment are still poorly understood. In a recent study in nonhuman primates (NHPs), we presented data suggesting that macrophages in different granuloma microenvironments have different polarization states that may control pathology through regulated expression of pro- and anti-inflammatory signals (5). We also showed that a characteristic spatial separation between pro- and anti-inflammatory macrophages could be observed in different granuloma types (5). Experimental studies have suggested that the presence of a dynamic balance between pro- and anti-inflammatory cytokines is necessary for achieving proper granuloma formation and control of M. tuberculosis infection (4, 56, 74, 75). In the context of macrophage activation, how this dynamic balance evolves over the course of infection remains to be elucidated and has been difficult to understand experimentally.
In this work, we build on our existing computational model of granuloma formation and function (GranSim), by integrating macrophage polarization data from our NHP model of granuloma formation, in order to test how cytokine receptor signaling influences macrophage polarization and bacterial containment (71, 76–78). We define a “macrophage polarization ratio” as a metric to understand how cytokine signaling translates into pro- or anti-inflammatory polarization and plasticity of single macrophages in a granuloma, which in turn modulates cellular functions, including antimicrobial activity and cytokine production (Fig. 1A). Gene expression profiles are triggered by simplified cytokine signaling models of STAT1, STAT3, and NFκB activation (Fig. 1B to D). We then define a “granuloma polarization ratio” by averaging macrophage polarization ratios over an entire granuloma, thus defining a metric of polarization at the tissue scale. While this metric is theoretical in nature, it predicts key polarization hypotheses, directly linked to immune function, with respect to contained granulomas which can be experimentally tested.
Schematic representation of macrophage polarization and gene expression dynamics captured in the agent-based model. (A) Simplified cartoon that represents the macrophage signaling cascade of the three pathways used in the study to capture M1 (proinflammatory) and M2 (anti-inflammatory) cues. Each signal was quantified and used to calculate a ratio that was then used to modulate several immune functions of the macrophage (chemokine and cytokine [TNF, IL-10] synthesis, bacterial killing). The modulation is implemented as a linear interpolation (the value of the ratio is normalized and multiplied by the parameter value regulating the correspondent immune function). (B, C, and D) Representation of how gene transcription dynamics is modeled through the three pathways (STAT1, STAT3, and NF-κB, respectively). For each of the pathways, we have a nuclear species state (denoted with N) and a response species state (denoted with R). The details are given in Materials and Methods. (E) Different gene transcription dynamic regimes (fast [green], intermediate [blue], and slow [red]) that are recapitulated by our in silico model, for each species. The two parameters whose values varied were the signal degradation time (or length [δ]) and the response/signal strength (β) of the response.
First, our computational model reproduces the spatial organization of the M1- and M2-like macrophages observed in many NHP granulomas (79). Second, we show that the granuloma polarization ratio can be used as a predictive metric of infection outcome as early as 2 months postinfection. Third, our in silico granuloma data also suggest that the length of time that the NF-κB signal is turned on (i.e., the signal activation interval [δ]) is a key regulatory mechanism, determining whether a granuloma contains infection with a minimal amount of inflammation or leads to uncontrolled bacterial growth. These studies provide clues to the importance of driving polarization ratios of macrophages within granulomas toward a containment phenotype and suggest how gene transcription dynamics might be manipulated pharmacologically and targeted for therapeutic control of infection toward this goal.
MATERIALS AND METHODS
Multiscale agent-based model.We recently published a next-generation multiscale computational model of lung granuloma formation and function during M. tuberculosis infection where we captured immunological processes over four different biological scales in the lung: the model operates at the molecular, intracellular, and cellular scales and provides readouts at the tissue scale (71). The main body of the model, GranSim, is based on the cellular scale and is populated with macrophages and T cells as distinct agents within an agent-based model setting (77, 78). Submodels operating within GranSim capture events at the molecular scale using ordinary differential-equation models (ODEs) that describe receptor-ligand dynamics for both IL-10 and TNF occurring within a single cell (71, 76). Diffusion of IL-10 and TNF in tissue is also tracked at the cellular scale (details can be found in references 69, 71, 76, 80, and 81). In the current study, we used TNF, IFN-γ, and IL-10 as drivers for macrophage polarization and as modulators of macrophage immune functions, such as cytokine synthesis and bactericidal activity. These assumptions are supported by observed biological data (8, 49, 54–56).
The model environment represents a 2-by-2-mm section of lung parenchyma tissue. Our two-dimensional square lattice comprises 100 rows of 100 microcompartments, with each microcompartment measuring 20 by 20 μm (approximately the diameter of a macrophage). The lattice is initialized with resident macrophages and vascular sources for recruitment of macrophages and T cells to the site of infection. Recruitment of cells into lung through vascular sources is based on our previous model studies and experimental data (71, 76–78, 81). Infection is initiated by placing a macrophage in the middle of the lattice and infecting it (through phagocytosis) with a single intracellular M. tuberculosis cell (82). We simulate infection over a time span of 200 days, with rules and interactions solved using a 10-min time step. Details regarding implementation of cell movement and molecule diffusion and degradation, as well as receptor-ligand dynamics, can be found in references 69, 71, 76, and 81. In order to make bacterial dynamics vary within known heterogeneous growth rates (83), rates are assigned by randomly choosing extracellular and intracellular M. tuberculosis doubling times from experimentally identified ranges defined a priori. Extracellular M. tuberculosis growth rates are sampled from a uniform distribution (defined by minimum and maximum values inferred from data; see Table 1) at the beginning of the simulation and assigned to each of the 100-by-100 microcompartments. A similar sampling scheme is used to assign intracellular growth rates to bacteria growing inside new macrophages entering the grid. Ranges for intracellular and extracellular bacterial growth rates are shown in Table 1, together with other parameters used in the model.
List of parameters varied in the uncertainty and sensitivity analysisa
Models of STAT1, STAT3, and NF-κB dynamics.To study the molecular mechanisms responsible for driving macrophage activation and polarization, we used a coarse-grained version of the NF-κB signaling module we previously implemented (80) and then modified our existing multiscale model to account for the intracellular signaling pathways mediated by STAT1, STAT3, and NF-κB. Although we use STAT1, STAT3, and NF-κB signaling as proxies for IFN-γ, IL-10, and TNF signaling, respectively, we understand that these signaling factors can be involved in other pathways, including IFN-α/β (STAT1), IL-6 (STAT3), and TLR4, T cell-, B cell-, or IL-1 receptor (NF-κB) signaling. Future work will be performed to identify the relationship between these factors and their contributions to signaling pathways in the granuloma and to add them to our computational model. We employed simplified versions of recently published models for STAT1, STAT3, and NF-κB pathways (88, 90, 94). We use two representative intracellular species for each signaling pathway: a phosphorylated nuclear species (denoted with N) and a response species (denoted with R) (Fig. 1B to D). The phosphorylated nuclear species (STAT1N, STAT3N, and NF-κBN) represent the phosphorylated forms of STAT1, STAT3, and NF-κB that exist in the nucleus following ligand stimulation. These nuclear phosphorylated species lead to downstream signaling (including transcription factors and mRNA), where they are referred to as the response species (STAT1R, STAT3R, and NF-κBR) (95, 96, 97). In this context, inhibition of signal activation for the response species represents the combined deactivation of these signal transducers in the nucleus, events mediated by dephosphorylation (STAT1 and STAT3) (97) and coupling with inhibitory subunits (NF-κB) and degradation of signaling molecules by proteolytic processes (96). The phosphorylated nuclear species have a maximum level (defined by STAT1MaxLevel, STAT3MaxLevel, and NF-κBMaxLevel) to capture saturation of ligand-induced stimulation.
Although the intracellular pathways are simplified by representing only a nuclear phosphorylated and a response species, we are able to recapitulate three gene expression regimes observed in many signaling pathways—fast expression, intermediate expression, and stable expression (98, 99). This is made possible via altering two model parameters for each pathway (see Fig. 1E), the signal activation interval (δ) and the signal strength (β). The signal activation interval for each pathway represents the length of time the signal is turned on. This interval is controlled by a degradation rate of the response species in the simplified signaling model (see Fig. 1B to D). The next step is to link these pathways to the existing multiscale model at the scale of the molecular species IL-10, TNF, and IFN-γ. To do this, we link models of IL-10 and TNF at the molecular scale (71, 76) to the new intracellular models of STAT3 and NF-κB, respectively (see Fig. 1), while STAT1 is similarly linked to IFN-γ. In the model, IFN-γ is secreted by T cells and its impact is restricted locally at the immunological synapse (100); thus, it is not necessary to explicitly model IFN-γ dynamics at the receptor scale (as we do for IL-10 and TNF). Thus, in the model, T cells produce IFN-γ as a proxy for activating STAT1N in neighboring macrophages with a given probability (PSTAT1).
There are other potential partners for STAT1 and STAT3 signaling activation, including type 1 interferons for STAT1 and IL-6 for STAT3 (20, 101). We are using TNF and IL-10 as activation signals since our current in silico model already tracks these two cytokines in a detailed way; we can include additional cytokines as more biological data supporting their roles in TB become available. As we have done previously, we use Poisson processes to describe the probability (PSTAT3 and PNF-κB) of ligand-induced activation of STAT3N and NF-κBN based on rate parameters (kSTAT3 or kNF-κB), threshold values (τSTAT3 or τNF-κB), and the respective concentrations of IL-10 bound to IL-10R and TNF bound to TNFR1 (71, 76, 80):
Macrophage polarization states.In the previous version of GranSim, we captured macrophages in resting, activated, infected, and chronically infected states (71, 76, 80). To explore macrophage polarization, we introduce new macrophage states: M0, M1, M2, and M1M2. M0 macrophages represent the initial resident alveolar macrophages and any unpolarized or nonactivated macrophages. M0 macrophages have no capabilities to secrete cytokines or chemokines and perform limited bactericidal functions (similar to the resting macrophage state in the previous version of GranSim). The states M1, M2, and M1M2 can all be considered subtypes of the activated macrophage state in the previous version of GranSim. Macrophages characterized as M1 have been polarized via stimulation of STAT1 or NF-κB or both. Macrophages characterized as M2 have been polarized via stimulation of STAT3. M1M2 macrophages capture the wide spectrum of macrophages between M1 and M2 (20), wherein these macrophages have been stimulated with a combination of STAT1, NF-κB, and STAT3.
Macrophage (RMP) and granuloma (RGP) polarization ratios.To determine the polarization of a macrophage toward a proinflammatory (i.e., M1) or anti-inflammatory (i.e., M2) phenotype, we define a “macrophage polarization ratio” (RMP) (Fig. 1A and B). RMP is a dynamic ratio comparing the amount of proinflammatory signals (STAT1R and NF-κBR) received to the amount of anti-inflammatory signals received (STAT3R) on an individual macrophage basis:
We normalize the RMP between 0 and 1 based on the minimum and maximum values of STAT1R, STAT3R, and NF-κBR dictated by the intracellular signaling model parameters (see Text S1 in the supplemental material [“macrophage polarization ratio—RMP calculation”]). Thus, macrophages that have a high value of RMP are polarized toward the M1 phenotype, while macrophages with a low value of RMP are polarized toward the M2 phenotype. We use a threshold of 1 for RMP as a theoretical value to label a macrophage as either M1 (RMP = >1) or M2 (RMP = <1). In order to use RMP (a cellular-level measure) at the tissue scale (i.e., the granuloma level), we collect macrophage polarization ratios for all macrophages within a single granuloma and average them. We call this composite average the “granuloma polarization ratio” (RGP) and use it as a metric at the granuloma scale. To make the averages consistent across different granuloma simulations, we account for macrophages only within a typical lesion size (∼1.5-mm diameter) and exclude from the average all the M0 (unpolarized/unactivated) macrophages.
Linking immune function with macrophage polarization using RMP.We use RMP, the macrophage polarization ratio, to link multiple immune functions (secretion of chemokines, secretion of two cytokines, and bactericidal ability) directly to macrophage polarization values (Fig. 1A). At the extremes of polarization, M1 macrophages secrete high levels of TNF and chemokines (CCL2, CCL5, and CXCL9/CXCL10/CXCL11) and very low levels of IL-10, while M2 macrophages secrete high levels of IL-10 and low levels of both TNF and chemokines (8, 102). Considering that a broad spectrum of macrophage polarization exists in vivo (20), we capture this by modeling secretion of TNF, IL-10, and chemokines as linear functions of log(RMP) (see Fig. 1A, bottom). Studies show that M1 macrophages express high levels of iNOS and low levels of arginase (Arg) whereas M2 macrophages express high levels of Arg and low levels of iNOS (5, 8). The two species compete for arginine as a substrate, which in the case of iNOS is used to produce antimicrobial RNI species and in the case of arginase is used to convert arginase to urea and l-ornithine, an upstream precursor of collagen (103, 104). Macrophages can simultaneously express iNOS and Arg enzymes (5, 103), suggesting that macrophage functional characteristics are a consequence of the enzyme abundance in that cell and not simply of enzyme presence or absence. Thus, we model macrophage bactericidal capabilities, condensing the iNOS and Arg pathways into a single response, as a linear function of the log of RMP (see Fig. 1A).
Arginase expression in macrophages as a consequence of IL-10 signaling was identified by few earlier studies (93, 105, 106). For example, in the study described in reference 105, the researchers used BALB/c mice infected with Leishmania and found that IL-10 signaled through STAT3 and unregulated IL-4 receptor on macrophages. This makes these macrophages more sensitive to IL-4 and thus sensitive to upregulation of arginase 1 in response to low IL-4 levels, reminiscent of the situation we see in the tuberculous granuloma. We do not model IL-4- or IL-4-mediated signaling, and we know that the number of IL-4-expressing T cells in granulomas is quite small, so IL-10 is used to represent an M2-polarizing factor.
Granuloma imaging.All procedures involving cynomolgus macaques (Macaca fascicularis) were performed in accordance with Institutional Animal Care and Use Committee protocols at the University of Pittsburgh. Macaques were infected with M. tuberculosis (Erdman strain), monitored for development of disease, and humanely euthanized as previously described (7); all tissues were from animals subjected to necropsy procedures for other studies. Granuloma-containing lung tissue fixed in 10% neutral buffered formalin and embedded in paraffin was deparaffinized and processed as previously described (5). Tissues were stained for phosphorylated STAT1 (pSTAT1) or pSTAT3. Immunohistochemical procedures involving pSTAT1 or pSTAT3 used Tris-buffered saline (TBS) as a diluent or wash buffer. Tissue sections were blocked in 1% fetal bovine serum (FBS)–TBS before being stained with rabbit anti-pSTAT1 (clone D3B7; Cell Signaling Technology, Danvers, MA) and mouse anti-pSTAT3 (clone M9C6; Cell Signaling Technology). Slides were washed with TBS before application of anti-mouse biotin (Vector Laboratories, Burlingame, CA) and Alexa Fluor 488-conjugated anti-rabbit (Life Technologies, Grand Island, NY) secondary antibodies. Biotinylated pSTAT3 was stained with Alexa Fluor 546-conjugated streptavidin (Life Technologies) followed by staining for CD163 (Thermo Fisher Scientific) with antibodies labeled with Alexa Fluor 647 via the use of Zenon labeling reagents (Life Technologies). Coverslips were applied to slides with Prolong Gold mounting medium containing DAPI (4′,6-diamidino-2-phenylindole) (Life Technologies), and images were acquired with a FluoView 1000 confocal microscope (Olympus, Center Valley, PA). For granulomas too large to be imaged in a single 200× field, multiple overlapping fields were acquired and the individual images assembled into a composite image using the Photomerge feature of Adobe Photoshop CS4 (Adobe Microsystems, San Jose, CA).
Model validation, UA, and SA.Most model parameter values used here are taken from references 69, 71, and 76 and are given in Table 1. We used semiquantitative studies from the literature on STAT1 (88, 89), STAT3 (90–93), and NF-κB (107) to specify ranges for parameters related to the macrophage polarization ratio introduced in this work (see Table 1 for details). Additionally, we relied on uncertainty analysis (UA) and sensitivity analysis (SA) techniques. These techniques can be used to efficiently explore the parameter space to inform baseline behaviors of the system (UA), as well as to quantify how parameter uncertainty impacts model outcomes. This allows efficient identification of critical model parameters that drive model behavior (SA; reviewed in reference 108). Here we used Latin hypercube sampling (LHS) for UA and partial-rank-correlation coefficients (PRCC) for SA. The LHS algorithm is a stratified Monte Carlo sampling method without replacement. It was used to generate 500 unique parameter sets, which were simulated in replication 3 times, and the averages of the outputs were used to calculate PRCC values. The details on the parameters that were varied in our LHS experiments are shown in Table 1.
Uncertainty analysis is also used to select in silico granulomas that recapitulate typical NHP granuloma images, as well as to generate statistics related to bacterial (CFU) numbers per granuloma over time (7, 109). The in silico granulomas are obtained without biasing signal activation intervals (δ) and signal strengths (β) toward one of the three pathways (STAT1, STAT3, and NF-κB; see Fig. 1E). Throughout the uncertainty analysis, the same range was used for each of the three signal activation interval parameters δ ([0.0001, 0.1] min−1; see Table 1), as well as for the three signal strength parameters β ([1, 25] min−1; see Table 1). These ranges allow us to recapitulate the three typical expression regimes (i.e., fast, intermediate, and sustained expression) exemplified in Fig. 1E, as observed in many signaling pathway systems (94, 98). We used a standard t test to compare the impacts of the three different pathways (i.e., STAT1, STAT3, and NF-κB) on granuloma development and maintenance.
Our comparison of in silico and ex vivo granuloma images is by no means a comprehensive comparison of all possible known histopathological outcomes seen in experimental studies in NHPs or in humans. However, we present a sampling here to provide evidence that our model can recapitulate a myriad of granulomas observed in NHPs, depending on the parameter values. We determined a baseline granuloma containment scenario (see Table 1 for baseline parameter values), where the CFU/granuloma grew in an almost uncontrolled manner in the first 4 to 6 weeks postinfection (as suggested by recent experimental data [109]) and then declined to lower bacterial levels, reaching a stable, controlled CFU level. We validated the model by performing virtual deletion simulations, where we set mechanisms (i.e., model parameters) to zero in upon infection to mimic experimental studies of gene knockouts (for example, IFN-γ and TNF; data not shown). All postoptimization analysis of data generated by our in silico model was done in Matlab.
Simulated granuloma classification.According to previously published NHP data (83, 109), at 4 weeks the median bacterial level in an individual granuloma is ∼2 × 104 CFU, with a large variability (between 1 × 102 and 1 × 106 CFU per lesion). Later time points show a decline of the average CFU level per granuloma (see references 83 and 109 for details), regardless of how the NHP has been classified with respect to infection outcome (i.e., latent or active TB). Although our in silico model tracks infection progression at only the single-granuloma level, we use bacterial load as a marker of granuloma outcome. If a granuloma by 200 days postinfection has a bacterial load below 2 × 104 CFU (e.g., if the CFU declines from a peak at day 30), it is classified as contained (containment scenario). If the bacterial load is above 2 × 104 CFU, then the granuloma is classified as disseminating (dissemination scenario).
RESULTS
Characteristic spatial distribution of polarized macrophages in TB granulomas.We first tested whether NHP data and model simulations gave similar spatial organizations with respect to macrophages and their phenotypes. Immunohistochemistry was performed on three granulomas from three different M. tuberculosis-infected NHPs. These granulomas were stained for pSTAT3 (pink) and pSTAT1 (green) (Fig. 2A, C, and E). The white areas at the center of the granulomas in panels A and C represent a necrotic core. The NHP granuloma images show a separation between pSTAT1-positive macrophages (green; localized closer to the center) and pSTAT3-positive macrophages (pink; localized at the outer regions of the granuloma). This supports previously published results from our group identifying the spatial organization of iNOS and arginase production by macrophages and showing that macrophages with lower iNOS expression, relative to arginase (M2-like) expression, localize to the granuloma's outer cuff whereas macrophages with higher iNOS (M1-like) expression were closer to the inner core (5). Representative in silico granulomas (chosen from our uncertainty analysis simulations) are shown in Fig. 2B, D, and F. Our model recapitulates this spatial organization, identifying granulomas that are necrotic (Fig. 2B and D) or nonnecrotic (Fig. 2F). As described in the Fig. 1 legend and in Materials and Methods, macrophages are classified based on their polarization ratio and labeled either M1 (green) if the polarization ratio is greater than 1 or M2 (pink) if it is smaller than 1. This spatial separation between M1- and M2-like macrophages is shared across different granuloma outcomes and is not due to a specific parameter set. Since at present we model only two cell types, e.g., macrophages and T cells, macrophages emerge as the largest fraction of cells that constitute our in silico granulomas. That is why a larger pink area is present in our in silico snapshots (Fig. 2B, D, and F) compared to the NHP granuloma images (Fig. 2A, C, and E), which also comprise other cell types (not labeled) such as neutrophils, B cells, and NK cells (4).
Comparison between granulomas from nonhuman primates and in silico granulomas generated by our model. (A, C, and E) Stains of 3 different granulomas for p-STAT3 (pink) and p-STAT1 (green). Green, M1-like macrophages; pink/purple, M2-like macrophages. (B, D, and F) Computer model snapshots at day 200 postinfection. As described in the Fig. 1 legend and in Materials and Methods, macrophages were classified as representing either an M1 (proinflammatory) or an M2 (anti-inflammatory) phenotype based on the polarization ratio. We label a macrophage as M1 (green) if its polarization ratio is greater than 1 and M2 (pink) if it is smaller than 1. Other cell types are as follows: effector lymphocytes (proinflammatory IFN-γ-producing T cells, light pink; cytotoxic T cell, purple; regulatory T cells, light blue), extracellular bacteria (olive green), vascular sources (gray), and necrotic spots (center; white x's). The three granuloma snapshots are generated by the in silico model with three different parameter settings (chosen from our uncertainty analysis simulations).
Granuloma polarization ratio dynamics are predictive of granuloma outcome.To better understand how the overall polarization state of macrophages in a granuloma relates to bacterial load, we examined granuloma polarization ratios (RGP) and CFU per granuloma over time (Fig. 3). Each point in Fig. 3 represents an in silico granuloma at a particular time point as generated in the set of virtual experiments described in Materials and Methods; each simulated granuloma is plotted over multiple time points. The containment and dissemination scenarios followed similar trajectories for bacterial growth in the first month postinfection (Fig. 3A and B), but while the containment scenarios showed a contraction in the levels of CFU/granuloma (Fig. 3C), the levels in the dissemination scenarios continued to increase over time (Fig. 3D). We then tested whether granulomas that contained infection had different granuloma polarization ratio dynamics with respect to the dissemination outcome. Figure 4A shows average trajectories of granuloma polarization ratios (RGP) for the two scenarios described in the Fig. 3 legend. The average time courses of RGP showed similar dynamics from the beginning of infection until day 60 (i.e., no significant differences between the two average trajectories), followed by a spike toward increased M1-like polarization ratios in the set of contained granulomas as early as 70 to 80 days postinfection. Later, the two trajectories crossed, and after day 170, contained granulomas showed significantly lower M1-like polarization. This suggests that a macrophage polarization ratio biased more strongly toward M1 phenotypes early on (i.e., within 2 to 4 months postinfection) appears necessary to contain bacterial loads at the single-granuloma level whereas at later time points (greater than 6 months) it may not be protective and may be detrimental due to excessive inflammation. Figure 4B shows the time course of TNF/IL-10 (calculated as the average value of ratios between TNF and IL-10 molecules) for the containment and dissemination scenarios (as shown in Fig. 3). Overall, the in silico model predicts that TNF levels are always higher (on average, 300 to 600 times higher) than IL-10 levels, whether the granuloma is contained or not. We also predict that average levels of TNF and IL-10 are 2 to 3 times higher in dissemination (data not shown). TNF/IL-10 ratios peak early during infection (first 4 weeks) and then decline rapidly within the first 6 weeks postinfection. After 6 weeks, the trajectories of the TNF/IL-10 ratios increase (the values corresponding to contained granulomas are always higher). However, by ∼10 weeks, while the contained granulomas slightly contract and stabilize their TNF/IL-10 levels, disseminating granulomas show sustained growth (mainly due to an increase in TNF; data not shown). Our model suggests that the higher TNF/IL-10 ratio in the first 40 to 80 days (Fig. 4B) drives a higher granuloma polarization ratio in the containment scenario (Fig. 4A).
CFU/granuloma versus granuloma polarization ratio (RGP), calculated as described in Materials and Methods. Each point represents an in silico granuloma at different time points during infection. The x axis shows the RGP versus its corresponding CFU value (y axis). Both axes are shown on a log scale. The data set shown is from an LHS of sample size 1,500 (n = 500 with 3 replications), with parameter ranges given in Table 1. Time is represented by labeling the data with different colors (see key in figure). The panels in the left column represent contained granulomas (panels A and C [labeled “containment”]). The panels in the right column represent granulomas that are not contained (panels B and D [labeled “dissemination”]). The classification distinguishing between those that are and those that not contained is based on the total bacterial load per granuloma, which falls below 2 × 104 after day 200 postinfection. Panels A and B track granulomas from day 10 to 40. Panels C and D track granulomas from day 50 to 200.
Granuloma polarization ratio (RGP) and TNF and IL-10 dynamics. (A) Granuloma polarization ratios (RGP) over time for both contained granulomas (n = 377) and disseminated in silico granulomas (n = 1,119). For ease of illustration, error bars are not shown. The time points between day 70 and day 130 as well as after day 170 show strong statistically significant differences between the two trajectories (*, P < 1e-2; **, P < 1e-3; ***, P < 1e-4). (B) Average trajectories of the ratios between the numbers of TNF and IL-10 molecules on the 2-mm-by-2-mm virtual environment between the contained and disseminated granuloma scenarios (as shown in Fig. 3).
Larger populations of polarized macrophages emerge in the dissemination scenario.Macrophage phenotypes were also differentially represented between the disseminating and containment granuloma scenarios (Fig. 5). Overall, the total number of macrophages in the dissemination cluster of granulomas was larger (Fig. 5A), likely due to more cellular influx as a result of higher bacterial loads, which could lead to higher levels of inflammation. The dissemination scenario showed a significantly larger population of M1-like polarized macrophages (∼8% versus ∼4%; Fig. 5D), as well as of infected macrophages (Fig. 5C), starting when the granuloma polarization ratios began to diverge between the two granuloma outcome groups, at about 2 to 3 months postinfection (see Fig. 4A). The larger numbers of cells in the M1 and infected macrophage populations in the dissemination group was mirrored by a smaller fraction of unpolarized macrophages, representing ∼95% of the macrophage population in the containment group (Fig. 5B). The larger populations of M2 and M1M2 macrophages in the dissemination group suggest that the lower granuloma polarization ratios seen in the 60-to-120-day window were driven by these macrophage phenotypes. In contrast, almost no M2 or M1M2 macrophages were observed in the containment group (data not shown).
Comparison of macrophage time courses for the in silico model. The y axis values represent either the cell count or the fraction (%) of cells over the total. For each panel, two trajectories are shown with standard error bars. The two curves represent average values over time for the two scenarios represented in Fig. 3. (A) Total number of macrophages. (B) Fraction of unpolarized macrophages (M0) over the total number of macrophages. (C) Fraction of infected macrophages (MI) over the total number of macrophages. (D) Fraction of M1 macrophages over the total number of macrophages.
Similarly, T cell counts were ∼2-fold to 3-fold higher in the dissemination class than in the containment granuloma group (Fig. 6A). There was no statistically significant difference between the dynamics of the fractions of the three different T cell phenotypes captured in the model (data not shown). However, while the ratio between IFN-γ-producing cells (Fig. 6B) and cytotoxic T cells (Fig. 6C) was approximately 1:1 for containment granulomas, T cell responses in disseminating granulomas were biased toward a higher (almost 2:1) ratio of cytotoxic T cells to IFN-γ-producing T cells.
Comparison of lymphocyte time courses for the in silico model. The y axis values represent the cell count. For each panel, two trajectories are shown with standard error bars. The two curves represent average values over time for the two scenarios represented in Fig. 3. (A) Total number of lymphocytes. (B) Total number of IFN-γ-producing lymphocytes (Tγ). (C) Total number of cytotoxic T lymphocytes (Tcyt). (D) Total number of regulatory T cells (Treg).
NF-κB signal activation dynamics characterize granuloma outcome.To identify the molecular mechanisms behind these predictions, we tested whether there is any significant difference in the dynamics of gene transcription parameters between the containment and dissemination granuloma groups. We first compared the signal activation intervals (δSTAT1, δSTAT3, δNF-κB) and signal activation strengths (βSTAT1, βSTAT3, βNF-κB) of the two granuloma groups shown in Fig. 3. In this context, shorter signal activation intervals can be interpreted as representing inhibition of the corresponding signal activation, since the length of time the signal is turned on (i.e., the signal activation interval) is regulated by proteolytic degradation or inhibition of the response species in the simplified signaling models (see Fig. 1B to D). We found significant differences between the signal activation intervals for the two different granuloma outcomes. Our in silico model predicts that granulomas that control bacterial loads have shorter NF-κB signal activation intervals than granulomas that are not able to control bacterial growth (Figure 7A [P = ∼1e-8]): δNF-κB is ∼50% stronger in the containment scenario than in the dissemination scenario (data not shown). We then compared signal activation intervals (δSTAT1, δSTAT3, δNF-κB) and signal activation strengths (βSTAT1, βSTAT3, βNF-κB) within the same granuloma groups. Figure 7B shows significant differences for the containment group. We predict that contained granulomas have NF-κB signal activation intervals that are significantly shorter than the STAT1 (δNF-κB = ≫δSTAT1 [P = ∼7e-6]) and STAT3 (δNF-κB = ≫δSTAT3 [P = ∼0.0115]) signal activation intervals, with the STAT3 intervals only marginally shorter than the STAT1 intervals (δSTAT3 = >δSTAT1 [P = ∼0.05]). Almost opposite dynamics were observed for the disseminating granuloma sets (Fig. 7C), with STAT1 signal activation intervals significantly shorter than the NF-κB intervals (δNF-κB = ≪δSTAT1 [P = ∼0.0018]) and STAT3 intervals only marginally shorter than the NF-κB intervals (δSTAT3 = >δNFκB [P = ∼0.05]).
Comparisons of gene transcription dynamics between and within the granuloma scenarios illustrated in Fig. 3. Only the signal activation interval (δ) parameters are shown, since they are the only significant ones. (A) Comparison of NF-κB signal activation intervals between containment and dissemination. The NF-κB signal activation interval is shorter (less stable) in containment than in dissemination (P < 1e-8). (B) Comparison between the signal activation intervals for the three different pathways in the containment scenario. The NF-κB signal activation interval is shorter (less stable) than the STAT1 interval (P < 8e-6) and the STAT3 interval (P < 0.012). (C) Comparison between signal activation intervals for the three different pathways in the dissemination scenario. The NF-κB signal activation interval is longer (more stable) than the STAT1 interval (P < 2e-3) and the STAT3 interval (P < 0.05).
Results of sensitivity analysis of the model also highlight the importance of NF-κB-related parameters as key mechanisms influencing most of the relevant readouts during infection progression and granuloma outcome (e.g., NF-κB signal activation interval, TNF-dependent NF-κB activation dynamics; see Table S1 in the supplemental material for details). In particular, NF-κB signal activation interval emerges as an important regulatory mechanism. Shorter NF-κB signal activation times play a key regulatory role both early and late during infection. The early role allows a prompt response to infection (shown by the spikes in the TNF/IL-10 ratio trajectories in Fig. 4B). Later during infection, a shorter NF-κB signal activation interval determines a faster contraction of the “activation” signal that is no longer needed, after the bacterial load is contained. Longer NF-κB signal activation intervals increase total TNF levels and the number of M1-like macrophages, inducing more inflammation and, consequently, more recruitment. This results in a large total number of macrophages at the site and larger lesion sizes. Overall, longer NF-κB signal activation intervals during the second and third month postinfection further sustain inflammation and cell recruitment, allowing the infection to persist and disseminate. TNF levels are higher for the dissemination trajectory late during infection, while they are similar for the two scenarios early on.
Interestingly, STAT1 signal activation interval lengths showed an impact almost opposite that of NF-κB on many of the readouts mentioned above, especially on the total TNF, total M1, and apoptosis levels (see Table S1 in the supplemental material for details).
These model results suggest that both NF-κB and STAT1 must be tightly regulated, with a shorter signal duration of NF-κB than of STAT1, thus preventing uncontrolled macrophage activation. Manipulating signal activation intervals can interfere with signaling dynamics and ultimately hamper and inhibit a protective response, as also shown in recent in vitro experimental and modeling studies (26, 27, 98, 99, 110). Our results show that the same principle can affect granuloma outcome in silico and, potentially, in vivo.
Spatial organization of polarized macrophages does not impact granuloma outcome.We investigated the possibility of a correlation between the spatial distribution of macrophages within a granuloma and granuloma outcomes. Figures S1, S2, and S3 in the supplemental material show the spatial distributions of macrophage polarization ratios over time between the two macrophage sets captured in Fig. 3. Each point represents a single macrophage, with its distance from the center of the granuloma on the x axis and its corresponding macrophage polarization ratio on the y axis. Excluded from these analyses are all unpolarized/unactivated macrophages (M0, as described in Materials and Methods). No clear differences emerged between the two outcomes within the first 100 days postinfection. After 3 to 4 months, a larger population of polarized macrophages was present in the dissemination scenario. However, this larger population did not show a spatial distribution different from that seen with the granulomas that control bacterial levels. In both scenarios, macrophages closer to the granuloma core always had a macrophage polarization ratio (i.e., they were near the site of infection and were getting the most M1-like stimulus) higher than that seen with the ones in the outer regions.
DISCUSSION
In pulmonary TB, granulomas in lungs are the key feature of infection. The composition of cells within primate granulomas is relatively well characterized (4), but their functions and the factors leading to different granuloma outcomes are not. Therefore, a better understanding of the cells and cytokines that both comprise and function within each granuloma will aid our understanding of both infection progression and outcome and provide possible targets for intervention to improve outcome. Here we used a computational-biology approach to further characterize the role of macrophage polarization within granulomas, described as M1, M2, or M1M2. Our modeling work builds on our recent experimental studies examining iNOS-producing (M1-like) and arginase-producing (M2-like) macrophages within granulomas from NHPs in which we found that many macrophages produced both enzymes but did so in differing amounts (5). In that work, we also characterized the spatial distribution of macrophages in granulomas, with epithelioid macrophages surrounding the necrotic core producing more iNOS than arginase (e.g., M1) whereas those in the outer cuff of the granuloma produced higher levels of arginase than iNOS (M2). The computational model used here allows us to perform focused in silico experiments that cannot be performed in vivo, thereby permitting exploration of potential immunomodulatory strategies for enhancing macrophage-based mycobacterial killing.
Our computational model replicates the spatial organization of M1-like and M2-like macrophages seen in many necrotic and nonnecrotic nonhuman primate (NHP) granulomas (see reference 5 and Fig. 2). However, this spatial organization does not appear to be correlated with granuloma outcome. This indicates that the spatial organization of polarized macrophages in NHP and in silico granulomas (M1-like macrophages closer to the core, with M2-like macrophages in outer granuloma regions) is a general property emerging from the host-pathogen interaction during granuloma formation that does not directly affect granuloma outcome. Recently published experimental and modeling studies (5, 71) indicate that a continuous and dynamic balance between pro- and anti-inflammatory signals is needed to control inflammation and restrain bacterial growth. Other modeling studies by Day et al. (111) and our group (69) defined the concept of “switching time” as the time necessary to switch from an M2-dominated (i.e., non- or anti-inflammatory) to an M1-dominated (i.e., proinflammatory) lung environment. In those studies, the authors elaborated on the biological relevance of increasing switching times in the context of M. tuberculosis infection, speculating that a delay in M1 polarization in the lung may be responsible for the bacterium gaining an initial “foothold.” Both studies described M. tuberculosis infection progression in the lung and predicted an optimal time window for a switch to occur in order to control bacteria dissemination. Here we tested potential mechanisms that drive macrophage polarization and plasticity and how the relative contributions of these broadly classified M1 and M2 programs evolve over time to initiate a protective immune response at the scale of a single granuloma.
To study how macrophage polarization in the granuloma evolves during infection and granuloma formation over many months, we defined and tracked two quantities: a macrophage polarization ratio (RMP) and a granuloma polarization ratio (RGP). We showed how, by following temporal dynamics of granuloma polarization ratios, we could predict granuloma outcomes. As early as 2 to 3 months postinfection, a spike in granuloma polarization ratio dynamics occurs for a set of granulomas that are able to contain infection. Granuloma polarization ratios for disseminating granulomas remain flat, without any significant oscillation toward either an M1-like or M2-like phenotype. M1 polarization timing and magnitude for an entire granuloma are significant factors in containing the infection. Our predictions are in line with those of Day et al. (111), where they predicted switching times of ∼50 days. Here we showed how stronger M1 polarization 2 to 3 months postinfection is necessary to drive a stronger protective immune response that is better at containing bacterial proliferation. In fact, the spike shown in Fig. 4A for the granuloma polarization ratios in the contained granuloma set can be interpreted as representing a more efficient way to contain higher bacterial loads during the first 4 to 6 weeks postinfection (83). We hypothesize that containment at this time point is critical because limiting bacterial replication with an early spike in M1 polarization limits excessive TNF expression and potentially pathological inflammation at later time points.
Differences in NF-κB gene transcription dynamics between the containment and dissemination granuloma scenarios represented in Fig. 3 suggest possible key regulatory mechanisms in driving a more or less protective response. We found that granulomas in the containment set have shorter NF-κB signal activation intervals than those granulomas in the disseminating set. Recently, Bai et al. published work suggesting that inhibition of NF-κB activation (i.e., shorter signal activation intervals) decreases survival of M. tuberculosis in human macrophages (112). They used an inhibitor of IκBa kinase to inhibit NF-κB activation and found a significant decrease in the number of viable intracellular mycobacteria recovered from THP-1 macrophages at both 4 and 8 days after infection. The granuloma environment is likely different from the environment described in these in vitro studies compared to the NHP or computational model; however, their approach highlights the pharmacological possibility of targeting an important key proinflammatory pathway in macrophages, such as NF-κB, and how this could be modulated in human infection to improve outcome. Another interpretation is that bacteria may actively interfere with NF-κB-mediated pathways to promote survival, as shown in other bacterial pathogen studies (36, 38, 39).
In our computational model, a way to achieve NF-κB inhibition is to decrease signal activation intervals (Fig. 1E). Increasing the length of these intervals translates into a more sustained signal, which can disrupt “normal” signaling dynamics driving an efficient inflammatory response; this has been demonstrated in both experimental and modeling studies of the NF-κB pathway (26, 98, 99, 110). During granuloma formation and development, there is a sustained level of TNF present. Thus, we speculate that controlling NF-κB activation (induced in the model by surface ligation of TNF by TNFR1) may be an important regulatory mechanism in vivo to contain infection, especially soon after an adaptive immune response arrives at the infection site and TNF-dependent M1 polarization is supported by IFN-y signaling (the two-signal concept that underpins classical macrophage activation). These results are in line with our previous study where we analyzed different NF-κB transcriptional responses, in the context of M. tuberculosis infection, including TNF secretion, chemokine secretion, inhibitors of apoptosis (IAPs), and macrophage-activating molecules (80). Here we show that with shorter NF-κB signal activation intervals, our in silico model produces stable necrotic granulomas, where most of the M. tuberculosis cells are slowly replicating or are nonreplicating, inducing less inflammation and lower levels of cell recruitment. Increasing NF-κB signal activation intervals results in extensive inflammation and, ultimately, uncontrolled CFU, as seen in the more complex pathologies observed in active TB (83). This might reconcile with the fact that some of the worst TB outcomes are associated with pneumonia or rapidly disseminating TB. Future work will integrate our higher-resolution molecular model of NF-κB (as shown in reference 80) into the current M1/M2 polarization implementation to get more mechanistic insights into the NF-κB pathway.
The hypotheses that we can test with our computational granuloma model, GranSim, have become increasingly complex as we have improved the model with additional laboratory-based data. This work addresses mechanistic issues that are otherwise difficult or impossible to address in vivo, including how the timing of macrophage polarization and polarization bias (toward M1 or M2) influences bacterial containment. Studying macrophages and T cell functions in granulomas is difficult in humans as it requires obtaining fresh lung tissue but is possible in cynomolgus macaques, the best animal model for recapitulating human TB (4, 113). Up to this point, most data have been acquired at an endpoint—removal of the granuloma as part of a necropsy. Newer technologies, including 18F-fluorodeoxyglucose (FDG)–positron emission tomography-computed tomography (PET/CT) imaging of humans (113, 114) and macaques (109, 115) with TB, has given us a new appreciation for the dynamic process of inflammation in the granuloma. Our in silico model complements these technologies by allowing us to test mechanistic hypotheses at the individual-granuloma level, including our hypothesis that distinct macrophage polarization and granuloma polarization trajectories are associated with differential control of bacterium and granuloma outcomes. Here we sought to investigate basic mechanistic factors underlying polarization of macrophage-based antimycobacterial responses using simplified models of STAT1, STAT3, and NF-κB signaling, through IFN-γ, IL-10, and TNF, respectively. Examining metrics such as macrophage and granuloma polarization ratios allows us to test whether macrophage-associated inflammatory factors are important and independent drivers of protective immunity and whether granuloma polarization ratios can describe the contribution of macrophage polarization to bacterial control in the granuloma. We anticipate that, as more biological data on macrophage function in the granuloma become available, we will be able to further refine the model, better understand interactions between host and pathogen, and develop new ways of shifting the therapeutic balance in favor of the host.
ACKNOWLEDGMENTS
We thank Paul Wolberg for the technical support during computational model development and Joe Waliga for website creation and management for the supplemental material for this paper.
Additionally, this research was supported in part by the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy's Office of Science. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under contract no. DE-AC02-05CH11231. This work also used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant no. ACI-1053575. This research was funded by NIH grants R01 EB012579 (D.E.K. and J.J.L. and J.L.F.) and R01 HL 110811 (D.E.K. and J.J.L. and J.L.F.) and R01 HL 106804 (D.E.K. and J.L.F.).
FOOTNOTES
- Received 14 August 2014.
- Returned for modification 10 September 2014.
- Accepted 28 October 2014.
- Accepted manuscript posted online 3 November 2014.
Supplemental material for this article may be found at http://dx.doi.org/10.1128/IAI.02494-14.
- Copyright © 2015, American Society for Microbiology. All Rights Reserved.
REFERENCES
- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.↵
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.↵
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.↵
- 84.
- 85.
- 86.
- 87.
- 88.↵
- 89.↵
- 90.↵
- 91.↵
- 92.↵
- 93.↵
- 94.↵
- 95.↵
- 96.↵
- 97.↵
- 98.↵
- 99.↵
- 100.↵
- 101.↵
- 102.↵
- 103.↵
- 104.↵
- 105.↵
- 106.↵
- 107.↵
- 108.↵
- 109.↵
- 110.↵
- 111.↵
- 112.↵
- 113.↵
- 114.↵
- 115.↵
