PLOS ONE RESEARCH ARTICLE A reaction-diffusion model to understand granulomas formation inside secondary lobule during tuberculosis infection Martı́ Català ID1,2, Clara Prats ID1,2*, Daniel López1, Pere-Joan Cardona2,3,4, Sergio Alonso ID1 1 Department of Physics, Universitat Politècnica de Catalunya, Barcelona, Catalonia, Spain, 2 Comparative Medicine and Bioimage Centre of Catalonia (CMCiB), Fundació Institut d’Investigació en Ciències de la Salut Germans Trias i Pujol. Badalona, Catalonia, Spain, 3 Experimental Tuberculosis Unit (UTE), Fundació Institut Germans Trias i Pujol (IGTP), Universitat Autònoma de Barcelona (UAB), Badalona, Catalonia, Spain, a1111111111 4 Centro de Investigación Biomédica en Red de Enfermedades Respiratorias (CIBERES), Madrid, Spain a1111111111 a1111111111 * clara.prats@upc.edu a1111111111 a1111111111 Abstract Mycobacterium tuberculosis (Mtb) is the causative agent for tuberculosis, the most extended infectious disease around the world. When Mtb enters inside the pulmonary OPEN ACCESS alveolus it is rapidly phagocytosed by the alveolar macrophage. Although this controls the Citation: Català M, Prats C, López D, Cardona P-J, majority of inhaled microorganisms, in this case, Mtb survives inside the macrophage and Alonso S (2020) A reaction-diffusion model to understand granulomas formation inside multiplies. A posterior chemokine and cytokine cascade generated by the irruption of mono- secondary lobule during tuberculosis infection. cytes, neutrophils and posteriorly, by T-cells, does not necessarily stop the growth of the PLoS ONE 15(9): e0239289. https://doi.org/ granuloma. Interestingly, the encapsulation process built by fibroblasts is able to surround 10.1371/journal.pone.0239289 the lesion and stop its growing. The success of this last process determines if the host Editor: Frederick Quinn, The University of Georgia, enters in an asymptomatic latent state or continues into a life-threatening and infective UNITED STATES active tuberculosis disease (TB). Understanding such dichotomic process is challenging, Received: May 12, 2020 and computational modeling can bring new ideas. Thus, we have modeled the different Accepted: September 3, 2020 stages of the infection, first in a single alveolus (a sac with a radius of 0.15 millimeters) and, Published: September 16, 2020 second, inside a secondary lobule (a compartment of the lungs of around 3 cm3). We have employed stochastic reaction-diffusion equations to model the interactions among the cells Copyright: © 2020 Català et al. This is an open access article distributed under the terms of the and the diffusive transport to neighboring alveolus. The whole set of equations have suc- Creative Commons Attribution License, which cessfully described the encapsulation process and determine that the size of the lesions permits unrestricted use, distribution, and depends on its position on the secondary lobule. We conclude that size and shape of the reproduction in any medium, provided the original author and source are credited. secondary lobule are the relevant variables to control the lesions, and, therefore, to avoid the evolution towards TB development. As lesions appear near to interlobular connective tis- Data Availability Statement: All relevant data are within the manuscript and its Supporting sue they are easily controlled and their growth is drastically stopped, in this sense secondary Information files. lobules with a more flattened shape could control better the lesion. Funding: CP, PJC and MC received funding from La Caixa Foundation (ID 100010434), under agreement LCF/PR/GN17/50300003; PJC received funding from Agència de Gesto d’Ajuts Introduction Universitaris i de Recerca AGAUR), Grup Unitat de Tuberculosi Experimental, 2017-SGR-500; CP, DL, Tuberculosis (TB) is an infectious disease that on 2017 killed 1.6 million people [1]. The same SA, MC received funding from Ministerio de year, nearly 10 million people developed the disease. Despite our best efforts, Mycobacterium PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 1 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Ciencia, Innovación y Universidades and FEDER, tuberculosis (Mtb) remains the bacteria able to cause the highest mortality by itself. Since with the project PGC2018-095456-B-I00. World Health Organization declaration of TB public health emergency in 1994 [2], TB death Competing interests: The authors have declared rate has been reduced from 23% in 2000 to 16% nowadays [1]. In 2015, the TB End Strategy that no competing interests exist. stated the goal of reducing its incidence by 50% in 2025 and 90% in 2035, as well as of reducing TB death to 75% in 2025 and 95% in 2035 [3]. It is estimated that between one quarter and one third of world population is infected with Mtb and that around a 10% of them will develop an active TB disease in the future years. Tuberculosis natural history has an extraordinary complexity and, in fact, there are still too many unknowns [4]. Tuberculosis infection starts when an Mtb is phagocyted by an alveolar macrophage (AM) at a pulmonary alveolus. Mtb resists bactericidal mechanisms induced by AM and replicates inside [5]. This capacity is the biggest challenge of Mtb and the main reason of its high mortality in humans along history. Under proper conditions, Mtb replicates approx- imately once a day [6]. When the intracellular bacterial load overcomes AM maximum tolera- bility, macrophage necrosis is triggered and thereby bacilli enter the extracellular milieu. These bacilli are phagocyted by other AM, which also fail to control the bacillary growth and are like- wise destroyed. This cycle is followed by a local inflammatory response and potentially ends once the specific immune response appears and controls it. This is the end of the progressive infection, which finally leaves an encapsulated TB lesion [6, 7]. Fibroblasts cells lead encapsula- tion with collagen, fibrin and other molecules. Fibroblasts cells are mainly located at pulmo- nary membranes like intralobular septae. In fact, this suggests that lungs structure may have an important role in TB infection dynamics, since initial infection distance to the nearest pulmo- nary membrane may determine the encapsulation capacity and, therefore, the final lesions size and bacillary load. Lungs are composed by around 2500 secondary lobules, and each lobule has an approximate volume of 3 cm3. There are two types of lesions: proliferative and exudative lesions. Proliferative lesions are mainly based in macrophages and lymphocytes, with intracellular bacilli inside the macro- phages and minimal necrosis, controled through time; exudative lesions are mainly based on necrotic tissue caused by neutrophilic infiltration where a lot of extracellular bacilli are accu- mulated [8–10]. If the response of the host is adequate, an Mtb infection can be completely cleared from the organism [11] or it can enter an asymptomatic latent state where the host is infected but not sick and cannot infect other people, corresponding to latent tuberculosis infection (LTBI). Nevertheless, if the immune and the inflammatory responses are not correctly balanced, the host can develop a more compromised situation, the active tuberculosis disease (ATB). Systems biology and computational models are great tools for increasing TB understanding [12]. Last years, several TB models have been built for a better understanding of different pro- cesses related with TB natural history [4]. In particular, GranSim [13] is a hybrid model that, under several modifications, is able to reproduce granulomas formation [14], encapsulation [15] and study drug resistance [16], among others. Bubble model [17, 18] is an Agent-based Model that reproduces the formation and growth of lesions, as well as the coalescence of neighboring lesions and their spreading by means of bronchial endogenous reinfection pro- cesses. Other mathematical models study different states in macrophage-bacilli competitivity [19]. In this work we develop a mathematical model to reproduce Mtb infection dynamics in an alveolus by means of Stochastic Differential Equations (SDEs). It is then generalized to a three- dimensional stochastic reaction-diffusion model that incorporates the space to reproduce Mtb infection, lesions formation and encapsulation in a secondary lobule. Finally, we identify the factors that determine lesions size and encapsulation time. PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 2 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Materials and methods Model structure Two different approaches are considered: the modeling of a single alveolus and the modeling of a secondary lobule. The difference between both models is the dimensionality of the space where the simulation takes place. Single alveolus system is a non dimensional model, therefore modeled by ordinary differential equations, and the secondary lobule is a three dimensional system modeled by partial differential equations. Next we explain the interactions among the different cells and elements of the model. We consider 5 types of cells in this model: bacilli (b), macrophages (m), neutrophils (n), T cells (T) and fibroblasts (f). Bacilli are distinguished between intracellular (bI), i.e. bacilli inside a mac- rophage, and extracellular (bE), i.e. bacilli in the extracellular milieu. Accordingly, macro- phages can be uninfected (mU), i.e. they do not have bacilli inside, infected (mI), i.e. they present bacilli inside but are not able to eliminate them, or activated-foamy (mAF), i.e. they present bacilli inside, are able to eliminate extracellular bacilli (activated) and drain intracellu- lar ones (foamy). Volume occupied by death elements (VO) and local inflammatory response (s) are also considered. Therefore, we consider the evolution of 8 types of cells and 2 additional processes. We take into account 16 processes, which lead to 10 reaction-diffusion equations, one for each element. These processes are described below. Note that the numbers in parenthe- ses relate to the description with the corresponding terms in the model equations and with the processes depicted in Fig 1. While the initial infection is closely related with the dynamics inside of an alveolus, the final part of the process, the encapsulation of the granuloma, takes into consideration the structure of the secondary lobule. For the first modeling appoach based on a single alvoulus, we study the effect of inflammatory and immune responses at the level of the alveolus, while for the sec- ond approach such responses are considered at a larger scale, the secondary lobule. For both approaches, we divide the whole evolution in four differential modules: Initial infection. Mtb infection starts when an AM phagocytes a bacillus in the alveolar surface (3, phagocytosis by non-infected macrophages). In this process, the extracellular bacilli become intracellular and non-infected macrophage becomes infected. Mtb resist bactericidal mechanisms induced by AM and replicates inside it (1, intracellular bacilli growth), as they use infected macrophages as sustenance. When intracellular bacillary load overcomes macrophage maximum tolerability (α), macrophage necrosis is triggered and thereby bacilli enter the extra- cellular milieu (5, macrophages’ lysis). Intracellular bacilli released become extracellular and the infected macrophage is eliminated. Then, other AM from closer alveoli enter the infected alveo- lus and the bacilli growth-macrophage lysis process occurs again. Infected macrophages also phagocyte extracellular bacilli (4, phagocytosis by infected macrophages) so that these extracellular bacilli also become intracellular. An explanatory schematic of this module is shown in Fig 1A. Inflammatory response entrance. Macrophages death and bacillary load’s increase trig- ger inflammatory response (7, inflammatory response activation). The variable s is increased from its initial value, 0, up to 1 after 3-5 days. This value is an approximation based on experi- mental observations [20]. When s = 1, inflammatory response is fully active. This increase in s leads to an entrance of neutrophils and uninfected macrophages in the alveolus (6, inflamma- tory entrance). The presence of neutrophils may be used by extracellular bacilli as a sustenance to grow (2, extracellular bacilli growth). There is a limit of how many extracellular bacilli can grow on a neutrophil (δ). If bacillary load decreases and infected macrophages disappear, inflammatory response is inhibited (8, inflammatory response inhibition). Macrophages and neutrophils cascade is unable to contain the infection, in fact, bacillary load grows exponen- tially. An explanatory scheme of this module is shown in Fig 1B. PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 3 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Fig 1. Four principal model’s modules. Interaction between different model elements. (A) Initial infection: interaction between bacilli and macrophages. (B) Inflammatory response: triggered by the presence of extracellular bacilli. (C) Immune response: infected macrophages trigger the mechanisms that activate this response. (D) Encapsulation: led by fibroblasts, it takes place when a lesion makes contact with a pulmonary membrane. Green arrows: activation; red arrows: inhibition; black arrows: transition; dotted arrows: chemo-taxis. Numbers identify specific processes: 1 intracellular bacilli growth, 2 extracellular bacilli growth, 3 phagocytosis by non-infected macrophages, 4 phagocytosis by infected macrophages, 5 macrophages’ lysis, 6 inflammatory response entrance, 7 inflammatory response activation, 8 inflammatory response inhibition, 11 immune response entrance, 12 macrophages activation, 13 extracellular bacilli elimination, 14 intracellular bacilli drainage and 15 fibroblasts chemo- taxis. See model’s details in Eqs (1)–(10). Part of this figure is adapted from [6]. https://doi.org/10.1371/journal.pone.0239289.g001 PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 4 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Immune response. Some infected macrophages trigger dendritic cells that travel to Lymph node to trigger immune response. Immune response consists of a T cells flow to the affected alveoli (11, immune response entrance). Infected macrophages are activated by contact with T cells (12, macrophages activation). This activated macrophages are able to phagocyte extracellular bacilli and eliminate them (13, extracellular bacilli elimination). Activated macro- phages have a short lifespan. When they are not able to eliminate more bacilli they convert to foamy macrophages [21–23]. These foamy macrophages leave the alveolus by draining the intracellular bacilli inside them (14, intracellular bacilli drainage). All the elements have a cer- tain lifespan; when they die, they leave a corpse that occupies a fraction of the space (9, death rate). A scheme of this module is shown in Fig 1C. Even when the lı́ving-time of the neutro- phils is limited (between 5 and 135 hours [24]), as they made neutrophilic extracellular nets (NETs) [25], and become part of the necrotic tissue, we have extended significantly the time of degradation. An explanatory scheme of this module is shown in Fig 1C. Encapsulation. During all these processes all elements diffuse, therefore occupying sur- rounding alveoli (16, diffusion). All elements diffuse by their own except intracellular bacilli, whose diffusion is conditioned by infected macrophages that contain them. The initial infec- tion forms a lesion due to cell accumulation. When lesions size is big enough to get in contact with pulmonary lobule membranes (interlobular connective tissue), encapsulation process is triggered. Encapsulation is led by fibroblasts cells that chemo-tact macrophages gradient (15, fibroblast chemo-taxis). Initially, fibroblasts cells are located at pulmonary lobules membranes. An explanatory scheme of this module is shown in Fig 1D. In Fig 1, a summary of the relations between the 10 variables of the model are shown, together with the different interactions involved, organized in the four modules previously mentioned. The corresponding 10 reaction-diffusion equations of the model are: � �3 ! � �3 ! d bI bI b ¼ mbI 1 þ ggmU bE þ ggmI bE 1 dt I amI |fflfflffl{zfflfflffl} amI |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ð3Þ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ð1Þ ð4Þ ð1Þ � � a b aðbbI þ ð1 baÞmI Þ lb bI m þ r gDm I rmI ; |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |{z} tm AF mI ð5Þ ð9Þ |fflffl {zfflffl } |fflfflfflfflfflfflfflfflfflfflfflfflffl {zfflfflfflfflfflfflfflfflfflfflfflfflffl } ð14Þ ð16Þ � �3 ! � �3 ! d bE bI b ¼ mbE 1 ggmU bE ggmI bE 1 dt E dn |fflfflffl{zfflfflffl} amI |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ð3Þ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ð2Þ ð4Þ ð2Þ þaðbbI þ ð1 baÞmI Þ lb nmAF g þ rðgDb rbE Þ ; |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |ffl{b zfflE} |fflffl{zfflffl} |fflfflfflfflfflfflffl {zfflfflfflfflfflfflffl } ð5Þ ð9Þ ð13Þ ð16Þ � � d b m ¼ ggmU bE þ k1 E sg lm mU þ rðgDm rmU Þ ; dt U |fflfflffl{zfflfflffl} m |fflffl{zfflffl} |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} ð3Þ ð3Þ |fflfflfflfflffl {zfflfflfflfflffl } ð9Þ ð16Þ ð6Þ d m ¼ ggmU bE ðbbI þ ð1 baÞmI Þ xmI Tg þ rðgDm rmI Þ ; ð4Þ dt I |fflfflffl{zffl fflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflffl {zfflffl } |fflfflfflfflfflfflfflffl {zfflfflfflfflfflfflfflffl } ð3Þ ð5Þ ð12Þ ð16Þ PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 5 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection d 1 m ¼ xm Tg m þ rðgD rm Þ ; dt AF |fflffl {Izfflffl } tm AF |fflfflfflfflfflfflfflfflfflm{zfflfflfflfflfflfflAF fflfflffl } ð5Þ ð12Þ |fflffl {zfflffl } ð16Þ ð14Þ � � d b 1 n ¼ k2 E sg n þ rðgDn rnÞ ; dt m tn |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} ð6Þ |fflfflfflfflffl {zfflfflfflfflffl } |{z} ð16Þ ð6Þ ð10Þ d T¼ lT T þ kmI sg þ rðgDT rTÞ ; ð7Þ dt |{z} |fflffl{zfflffl} |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} ð9Þ ð11Þ ð16Þ d f ¼ rðgwf jrgjrf Þ ; dt |fflfflfflfflfflfflfflfflfflffl {zfflfflfflfflfflfflfflfflfflffl } ð8Þ ð15Þ d 1 VO ¼ ðbbI þ ð1 baÞmI ÞVm þ lm mU Vm þ lT TVT þ nVn ; dt |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflffl {zfflfflfflffl } |fflffl {zfflffl } tn ð9Þ ð5Þ ð9Þ ð9Þ |fflffl{zfflffl} ð10Þ d s s ¼ rU mI bE ð1 sÞ rD ; dt |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} bE þ 1 ð10Þ ð7Þ |fflfflfflfflffl{zfflfflfflfflffl} ð8Þ where g is the unoccupied fraction of volume in each alveolus: ðmU þ mI þ mAF ÞVm þ nVn þ TVT þ VO g¼1 ; ð11Þ Va and k1(x) and k2(x) are the entrance rate of macrophages and neutrophils, respectively. They are quadratic functions that depend on extracellular bacilli—macrophages ratio. When this ratio is high enough, immune response is neutrophils-based (NBR) and most of the entering elements are neutrophils, otherwise, if macrophages are able to control extracellular bacilli, the immune response is macrophages-based (MBR) and most of the entering elements are macro- phages. The threshold that differentiates both is half of macrophages tolerability (α/2). The total amount of external agents flux (Y) is conserved: k1 ðxÞ � Vm þ k2 ðxÞ � Vn ¼ U; ð12Þ U 4�þ x2 þ � a2 k1 ðxÞ ¼ ; ð13Þ Vm 4x2 þ a2 U 4ð1 �þ Þx2 þ ð1 � Þa2 k2 ðxÞ ¼ : ð14Þ Vn 4x2 þ a2 Macrophages and neutrophils entrance rate profile per day can be observed in Fig 2 for values of extracellular bacilli and total amount of macrophages ratio between 0 and 2α. PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 6 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Fig 2. Inflammatory response profile. Profile of the inflammatory response as a function of extracellular bacilli and total amount of macrophages. (A) Number of cells that enter per day at the alveoli when inlfammatory response is fully activated. (B) Volume fraction of macrophages (blue) and neutrophils (green) that enter in alveoli due to inflammatory response. https://doi.org/10.1371/journal.pone.0239289.g002 Note that the 3D diffusion equations part is marked in blue. If the blue parts of the equation are not considered, we recover the equations for a single alveolus model. Parameters This model depends on a total of 27 parameters that are obtained from bibliography or adjusted from experimental observations. All the parameters, their baseline value and corre- sponding references are shown in Table 1. Parameters Y, κ and ξ are adjusted to reproduce Bru & Cardona observations [26]. Parameter δ is adjusted from Marzo et al figures [20]. Parameter ρ is adjusted from Marzo et al observations [20]. Parameters β and γ are computed using simple Agent-Based Models (ABM). Adjusting β. The parameter β is fitted by means of an ad-hoc ABM designed for studying the distribution of bI among mI and the resulting amount of macrophages’ lyses. We develop a simple ABM where macrophages are the agents with only one property: the number of bacilli inside them. All other elements behave like the model described before. There is no diffusion, as just one alveolus is considered (single alveolus model). Macrophage lysis is triggered when its bacillary load is higher than macrophages maximum tolerability (α). Macrophages death rate due to bacilli growth is computed at each time step. As a result of the simulation, a func- tion that relates the number of intracellular bacilli, number of infected macrophages and lysis rate due to intracellular bacilli growth can be fitted: lysisðmI ; bI Þ ¼ b1 bI þ b2 mI ; ð15Þ where βi are the adjusted constants. Fitting goodness R2 is computed and values over 0.8 are obtained for different parameters sets exploring diverse initial conditions. Using baseline val- ues and a set of random initial conditions we obtain: β1 = 0.01816 day-1, β2 = −0.1007 day-1, R2 = 0.895. Note that, as expected, the number of bacilli needed to kill one macrophage is α: b1 � a þ b2 � 1 � 1: ð16Þ PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 7 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Table 1. Parameters description and values. Parameter with source ABM is adjusted from Agent-Based Model, details are on the text. Parameter with source ABM� is fit- ted from ABM and corrections are implemented. Parameters with adj. source were adjusted to reproduce tuberculosis infection dynamics according to that source. Param- eters Db, Dm, Dn, DT and χf are marked with ? because they are not directly explored. These parameters are explored using D and DF parameters, related with Eq 24. Parameter D is explored between half and double its baseline value (4.6�10-9) and DF is explored between 0.5 and 4. parameter description values baseline value units explored values source μ macrophages growing rate 0.69–1.04 0.693 day-1 0.35-1.39 [6, 27] α macrophages maximum tolerability 32–64 60 - 30-120 [6] β bacilli kill macrophages 0.01–0.04 0.0182 day -1 adj. ABM γ phagocytation rate 8 �10-4–8 0.05 day-1 0.025-0.1 ABM� δ neutrophils capacity 10-40 30 - 15-60 adj. [20] κ T cells recuitment 0.8 day-1 0.4-1.6 [19, 26] ξ macrophages activation 0.01 day-1 0.005-0.02 [19, 26] ν activated killing rate 150 day-1 75-300 adj. [6] -4 ρU inflammatory activation 10 day-1 10-5-10-3 adj. [20] ρD inflammatory inhibition 1 day-1 no adj. [20] 5 3 -1 5 6 Y inflammatory flow 5 � 10 μm day 2 � 10 -10 adj. [6] �+ macrophages proportion in NBR 5/8 - no adj. [6] �− macrophages proportion in MBR 50/53 - no [28] λb bacilli death rate 1- 8 � 10-3 10-2 day-1 0.005-0.02 [29] λm macrophages death rate 0.023 day-1 0.01-0.05 [30] τm activated macrophages lifespan 3–10 5 day 2.5-10 [6, 31] tn neutrophils degradation time 100 day 50-200 adj. λT T cells death rate 0.20–0.33 0.33 day-1 0.16-0.66 [6] 3 Vm macrophage volume 5000 μm no [32] Vn neutrophil volume 300 μm3 no [33] VT T cell volume 200 μm3 no [34] Va alveolus volume 4 �106–6 �107 107 μm3 no [35] Db bacilli diffusion Dm/10 4.6 �10-10 cm2s-1 ? [31] Dm macrophages diffusion 10-5–10-11 4.6 �10-9 cm2s-1 ? [31] -5 -11 -9 Dn neutrophils diffusion 10 –10 4.6 �10 cm2s-1 ? [31] DT T cells diffusion 10-5–10-11 4.6 �10-9 cm2s-1 ? [31] χf fibroblasts chemotaxis coefficient 9.2 �10 -9 cm s 2 -1 ? adj. [6] https://doi.org/10.1371/journal.pone.0239289.t001 Then, imposing Eq (16) and redefining β1 as β we obtain a new expression for the lysis func- tion: lysisðmI ; bI Þ ¼ bbI þ ð1 baÞmI : ð17Þ Using baseline values we obtain β = 0.0182 day−1 and R2 = 0.884. The value of β is calculated at the start of every simulation and depend on the parameters set used. Parameters α and μ are the parameters that mostly determine β value. Adjusting γ. The parameter γ is also estimated using another ad-hoc ABM that simulates a simple macrophages-bacilli interaction in an alveolus. We consider a 2D system where bacilli are considered fixed rectangles of width 0.6 μm and length 4 μm [36] and macrophages are mobile circles with radius rm = 10.6 μm [32]. Simulations take place in a L × L square with periodic boundary conditions which represent alveolar surface, L = 532 μm [35]. Macrophages diffuse moving at constant velocity but with a changing direction. Using Fürth formula [37] a PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 8 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection diffusive movement can be characterized by the following equations: d ~ x ðtÞ ¼ ~ vðtÞ; ð18Þ dt ~ vðtÞ ¼ v0 � ð cos yðtÞ; sin yðtÞÞ; ð19Þ d yðtÞ ¼ s � ZðtÞ; ð20Þ dt where v0 is a constant velocity that depends on diffusive coefficient, Dm, and persistence time, τp, that for macrophages can be approximated as τp = 2000 s [38]: sffiffiffiffiffiffiffiffi 2Dm v0 ¼ ; ð21Þ tp and η(t) is white noise of amplitude σ that depends on persistence and time step, dt. sffiffiffiffiffiffiffiffi 2 s¼ : ð22Þ tp dt At each time step, macrophages move and phagocyte bacilli at a distance smaller than its radius, (rm). These bacilli are eliminated from the system. Every time that a macrophage phagocytes a bacillus, macrophage is stopped during tphag = 10 minutes [27]. At each time step, the number of bacilli, macrophages and phagocytations that take place are computed. Different random initial conditions are explored. As a result, the phagocyte ratio can be approximated as: phagocyteðm; bE Þ ¼ gmbE : ð23Þ Different sets of parameters adjust to different γ values. R2 is computed being always higher than 0.90. Using baseline values γ = 0.56 day-1 and R2 = 0.9882. This simplified model give us a range of values for γ, γmin = 8 � 10-4 day-1 (Dm = 10-11 cm2 s-1) and γmax = 7.6 day-1 (Dm = 10-5 cm2 s-1). Due to biological factors as alveolus saturation, long range interactions and occupa- tion factor we considered γ smaller and its baseline value was set to γ = 0.05 day-1. Implementation The model was implemented in MATLAB using Poisson τ-leap method to integrate differen- tial equations [39], τ = 0.1 days for the one alveolus model and τ = 0.001 days for the secondary lobule model. Diffusion was computed with finite differences, Δx = 0.3 mm. Transition proba- bilities were determined by the differential equation and the success or not was determined using Poisson random numbers to assure that all quantities are integers. Poisson random numbers are also used in finite differences. For the secondary lobule model (3D simulations) we consider a cubic secondary lobule of 3.4 cm3 (1.5 × 1.5 × 1.5 cm). Secondary lobule is formed by 125000 (50 × 50 × 50) alveoli, cubes of 0.3 mm side [35]. Reaction equations take place in each alveoli and diffusion is con- sidered with surrounding alveoli. The cube is surrounded by a thin membrane (1 alveolus thickness) where fibroblasts concentration is maximum. Membrane fibroblasts are considered inexhaustible (boundaries are always full of fibroblasts). PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 9 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Results Model of a single alveolus Two final states are observed: Proliferative and exudative scenarios. Single alveolus model results are computed eliminating diffusive (blue) parts of Eqs (1)–(10). Exploring differ- ent set of parameters and diverse initial conditions we see that there are two possible final states as it can be seen in Fig 3. (1) Proliferative scenario: bacilli population is eliminated by a flux of macrophages that are activated by T cells. Temporally, part of the alveolus is filled by macrophages, however, finally the number of macrophages returns to 1. Some of the space remains occupied by dead cells. (2) Exudative scenario: bacilli persevere and are able to kill all macrophages that try to contain the infection. Most of the present cells in this scenario are neutrophils thus favouring extracellular growth of bacilli. At the end, all bacilli are extracellular due to macrophages lysis and alevolus is occupied by dead cells of macrophages, neutrophils and T cells. Under certain values of the parameters, bistability of both states is observed and the fluctua- tions determine the final state. According to these different simulation possibilities we can define three zones of parameters: Proliferative zone where all simulations end in a proliferative scenario, bistable region where a simulation can end in exudative or proliferative scenario depending on the particular realization, and exudative zone where all simulations end in exu- dative scenario. For each set of parameters we can compute the probability to end in exudative state, pexu, or proliferative scenario, ppro. Note that pexu = 1 − ppro. Activation rate of macrophages determines transition between proliferative and exuda- tive scenarios. In Fig 3, a sample of simulations with three different set of parameters are shown, one set for each defined zone (proliferative, bistable and exudative). The parameter used to explore different parameters zone is macrophage activation rate (ξ). Different values used were: proliferative (ξ = 0.5 day-1), bistable (ξ = 0.15 day-1) and exudative (ξ = 0.01 day-1). For the set of parameters of the bistable zone, the probability to end in exudative scenario is pexu = 0.22. This probability is computed using the results of 10000 simulations till day 200. Initial behaviour is very similar for all sets of parameters, inflammatory response is trig- gered by first macrophages lysis that takes places 6.4 ± 2.3 days after initial infection and its fully awaken around day 16th. Firsts T cells arrive at the alveolus at day 14.2 (1 T cell), and around day 17th there are around 25. At day 16.7 ± 2.8, the first activated macrophage is observed, after 19th day more than 10 macrophages are active. Transition to exudative scenario depends on bacilli growth and T cells and macrophage activity. We perform a sensitivity analysis to determine the input parameters effect in out- come variables. We employ the methodology described in [40]. Parameters are sampled using a Latin Hypercube Sample technique [41] between the values specified in Table 1. We explore a total of 14 parameters, which are those that do not have a well defined value. In most of the cases, the exploration ranges are between half and double of their baseline value, but in some parameters of higher interest exploration ranges are extended. We use 1500 points to perform the exploration. Partial Rank Correlation Coefficient (PRCC) is computed between input parameters and outcome variables along time. In Fig 4 we show the correlation between six parameters (μ, α, δ, γ, Y and ρ) and the outcomes intracellular bacilli (bI), extracellular bacilli (bE), and total macrophages (m). Parameters κ, ν, ξ, λb, λm, λT, τm and tn are not shown in Fig 4 because they do not present correlations bigger than 0.2 with these outcome variables. It can be seen that as expected parameter μ has a positive correlation with the number of bacilli, long term correlations decrease because alveolus is filled. Initially μ reduces macrophages quantity due to bacilli growing inside, but later this effect triggers inflammatory response and macro- phages numbers increase. Parameter α contributes to increase bacilli quantities because more PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 10 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Fig 3. Single alveolus model simulations with three different set of parameters. (A, D, G, J, M) Proliferative set, (B, E, H, K, N) bistable set and (C, F, I, L, O) exudative set. For each set of parameters there can be seen the total amount of intracellular bacilli, in red (A, B, C), the total amount of extracellular bacilli, in orange (D, E, F), the total volume occupied by macrophages, in blue (G, H, I), the total volume occupied by neutrophils, in green (J, K, L) and the occupied alveolus volume fraction, in black (M, N, O). Thin and bright lines represent single simulations, 10 single simulations are shown as example for each set. Thick and dark lines are the mean value of 10000 simulations. Simulations were 200 days long, but only the first 140 are shown. https://doi.org/10.1371/journal.pone.0239289.g003 PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 11 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Fig 4. Sensitivity analysis for intracellular bacilli, extracellular bacilli and total macrophages. Partial Rank Correlation Coefficient (PRCC) is computed between a total of six input parameters (see legend) and outcome variables. A total of 1500 parameters sets are used for exploring the parameter space defined in Table 1 using Latin Hypercube Sample technique. A total of 14 parameters are explored, the ones that are not shown it is because their PRCC values is smaller than 0.2. Gray horizontal lines are marked at ±0.043 to mark significance threshold for PRCC values using a significance level of 0.05. PRCC are computed in each time step (dt = 0.1 day). For each set of parameters a total of 10000 simulations are done. https://doi.org/10.1371/journal.pone.0239289.g004 bacilli can duplicate inside a macrophage, these also trigger inflammatory response and more macrophages enter to the alveoli. Parameter δ increases the number of extracellular bacilli because more bacilli can duplicate over a neutrophil. Parameters Y and ρ have a very similar effect because both are related with inflammatory response (Y increases the entering flux and ρ activates faster the inflammatory response). This accelerates the process and more bacilli and macrophages are seen before. Parameter γ is important to determine how fast are bacilli engulfed and its effect is mainly seen after first macrophage lysis. As seen in Fig 3, the set of parameters determines final state probabilities. We compute PRCC between 14 input parameters used in sensitivity analysis and the probability to end in exudative scenario. In Fig 5 PRCC final values are shown. Increasing T cells death rate (λT) and macrophages growing rate (μ) or reducing macrophages activation (ξ), T cells recuitment (κ), activated macrophages killing rate (ν) and activated macrophages lifespan (τm) increase the probability to end in exudative scenario. Therefore, parameter ξ is not the unique parame- ter that determines the zone where a desired set of parameters belong. Model of secondary lobule Fibroblast encapsulation stops the expansion of the lesion in the model. Initially all alveoli present one uninfected macrophages except a certain alveolus that has an infected mac- rophage with a bacillus inside. The resulting evolution of the infection can be observed in Fig 6. We show the total bacilli amount (b), total macrophages amount (m) and fibroblasts con- centration (f). In supplementary material there is a video with lesions growth and encapsula- tion process. Initial infection starts with an infected macrophage and an intracellular bacillus that dupli- cates. Intracellular bacilli duplicate until lysis of first macrophage; then, approximately α extra- cellular bacilli are released and diffuse to nearby alveoli. Infection continues as seen in Fig 3 (exudative set) but extending to surrounding alveoli. It forms an spherical lesion with a radius that grows exponentially. This unstoppable process continues until lesion makes contact with secondary lobule membrane (septae), which occurs around day 51 as seen in Fig 6, where base- line values are used. This contact triggers the encapsulation process. Fibroblasts contain lesions PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 12 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Fig 5. Partial Rank Correlation Coefficient (PRCC) between input parameters and exudative probability. PRCC between ξ, κ, ν, λb, γ, tn, Y, λm, δ, α, τm and λT and exudative probability are shown in blue bars. Infection probability is computed at time = 200 day as the presence or not of bacilli (bI + bE). Gray horizontal lines are marked at ±0.043 to mark significance threshold for PRCC values using a significance level of 0.05. A total of 1500 parameters sets are used to explore the parameter space defined in Table 1 using Latin Hypercube Sample technique. For each set of parameters a total of 10000 simulations are done. Error bars are computed as the standard error of Spearman’s correlation coefficient. https://doi.org/10.1371/journal.pone.0239289.g005 growth and an encapsulated lesion is formed (day 73). Encapsulation process finishes around day 126 and no diffusion is observed after this time. The shape of the lesion is spherical because encapsulation process is faster than lesion growth. If encapsulation is slower (it can be seen with a different set of parameters), the shape of the lesion is oval. In fact, if fibroblasts chemo- taxis is slow enough, lesions may occupy the whole secondary lobule and encapsulation is not formed. Encapsulation time and volume of the lesions depends on fibroblast spreading. An uncertainty analysis is performed following the same procedure as in the single alveolus model. A total of 100 simulations are run, exploring 16 input parameters (Table 1). Diffusion coefficients are also varied using two parameters, D and DF. Values of diffusion coefficients are related following: wf Dm ¼ Dn ¼ DT ¼ 10 � Db ¼ ¼ D: ð24Þ DF Parameter D is varied between the double and its half of its baseline value and DF is explored between 0.5 and 4 times its baseline value. PRCC are computed between input parameters and outcome variables. Total amount of intracellular and extracellular bacilli, and macrophages are similar than the values seen in the single alveolus model, see Fig 4. In Fig 7A PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 13 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Fig 6. Evolution of 3D model. Evolution of the secondary lobule model in a cubic 50 × 50 × 50 alveoli grid. The first row is the 3D representation of the secondary lobule. Fibroblasts are shown in green and the occupied space (macrophages, T cells, dead cells, bacilli. . .) is in blue. Second, third and fourth rows represents a slice (z plane at the height of initial infection). Concentration of fibroblasts in shown in green (second row), the total number of macrophages (uninfected, infected and activated) in logarithmic scale in blue (third row) and, finally, the total number of bacilli (intracellular and extracellular) in logarithmic scale in red (fourth row). These quantities and representations can be observed at different days after initial infection, each column represents a different time of simulation. https://doi.org/10.1371/journal.pone.0239289.g006 we observe PRCC values between volume of the lesions and the explored parameters. Lesions volume is determined by two initial parameters D, diffusion of the elements, related with lesions growth, and DF, fibroblasts diffusion, related with encapsulation velocity. If the lesions grow faster it is more difficult to encapsulate them and final lesions are bigger. In fact, the faster the encapsulation speed, the smaller the lesion. This relation is observed in Fig 7B. Size of the lesions is also determined (less influence than D and DF) by α, μ and ρ. These 3 parame- ters increase bacillary load, as shown in Fig 4. This causes a higher diffusion of bacilli and mac- rophages, thus a bigger lesion. In Fig 7B there is PRCC between encapsulation time and input parameters. Encapsulation time is defined as the elapsed time between first time that fibroblast diffusion is non zero and the time where it becomes zero again. Parameter DF reduces encapsulation time because it increases fibroblasts diffusion and, as seen in Fig 7A, parameter D increases final lesions volume. We have observed that bigger lesions have a larger encapsulation time, then we would expect that D increases encapsulation time. However, in Fig 7B it is seen that, contra-intui- tively, it decreases. This is due to the fact that increasing D also increases fibroblasts diffusion but, in fact, as the mean value of DF is higher than 1 (it is 2.25) the increase in fibroblasts diffu- sion is amplified. Parameters DF and D are not the only parameters that determine lesions final volume. PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 14 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Fig 7. Sensitivity analysis of 3D model. (A) Partial Rank Correlation Coefficient (PRCC) values between volume of the final lesion and explored input parameters of the model. 100 sets of input parameters are used. Gray horizontal lines are marked at ±0.166 to mark the significance threshold for PRCC values using a significance level of 0.05. Error bars are computed as the standard error of Spearman’s correlation coefficient. (B) PRCC values between encapsulation time and input parameters. 100 sets of input parameters are used. Gray horizontal lines are marked at ±0.165 to mark the significance threshold for PRCC values using a significance level of 0.05. Error bars are computed as the standard error of Spearman’s correlation coefficient. (C) Distance to the nearest pleura against the volume of the resulting lesion. Each point is the mean value and standard deviation for three different simulations. Different values of D are explored. In black, baseline value, in red, baseline value divided by 2 and, in blue, baseline value multiplied by 2. (D) Volume, in blue, and radius, in orange, of the lesion using baseline values. Dotted black line is generalized logistic second degree curve approximation defined by Eqs 25–27. https://doi.org/10.1371/journal.pone.0239289.g007 Note that Fig 7 error bars are bigger than the ones observed previously in Fig 5, due to the low number of realizations done in the secondary lobule model, because it is more computa- tionally demanding. Distance to the nearest pulmonary membrane strongly determines final volume of the lesion. Different simulations are carried out with different initial infection positions. As cubic secondary lobules are considered, the distance to membranes can be computed from three different axis (x, y and z). The initial infection is considered in the middle of the z axis. Then, the considered distances are 7.5 mm for each one of the sides. The distance to y mem- brane is fixed at 3 mm for the closest side and 12 mm for the other side. The distance to x membrane is varied for the different simulations. The considered distances are: 0.6 mm, 1.5 mm, 2.1 mm, 3.0 mm, 4.5 mm and 6.0 mm. For each different distance, three simulations are run. All simulations are repeated considering a doubled baseline value for D and using the half of its value. In Fig 7C the volume of the lesion as a function of the distance to the nearest x membrane for the three different values of D can be seen. Two different behaviours can be distinguished. When the distance is smaller than 3 mm, a clear proportionality between distance and volume is observed. The bigger the distance the PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 15 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection bigger the volume. As small data is available, it can not be distinguished between exponential growth or a polynomial growth of a certain power. It is expected that volume is proportional to the cube of the distance to the nearest membrane. When the distance is bigger than 3 mm, this behaviour changes. This is due to that the near- est distance thought x axis (that is the plotted one) is not the small one because the nearest y membrane is at 3 mm. It can be seen that the volume does not increase now with the distance. Then it can be seen a clear dependence between the volume of the lesion and the distance to the closest membrane. In Fig 7C it is seen that the volume of the lesions increase with the diffusion (bigger values of D parameter). In blue simulations with the double of D value, in black the baseline value and in red half of the baseline value. All the colors follow the same shape. This result is also seen in Fig 7A. Radius of the lesion can be fitted with a logistic equation. In Fig 7D lesions growth is observed (volume in blue and radius in orange) for a single secondary lobule simulation using as initial infection point an alveolus that is at 2.7 mm of the closest membrane. This growth can be approximated by a generalized logistic function [42] of second degree for the radius: � �2 ! d rðtÞ rðtÞ ¼ v � rðtÞ � 1 : ð25Þ dt rmax That can be solved imposing that there is a minimum radius (rmin) at a desired time (tmin): r rmin rðtÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffimax ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffi : ð26Þ rmin þ ðrmax rmin 2 2 2 Þ � expð 2vðt tmin ÞÞ Catala et al [18] adjusted an agent-based model based on experiment observations of mini- pigs. A logistic growth was adjusted using the following parameters: rmax = 2.62 mm, v = 0.1324 day-1, rmin = 0.0207 mm and tmin = 14 days. Imposing tmin = 14 days to obtain com- parable results we can adjust the logistic as: rmax ¼ 2:62 mm; v ¼ 0:1324 day 1 ; ð27Þ rmin ¼ 0:0207 mm: Comparing both models it can be seen that in this model (human) lesions are assumed to be bigger but to grow slower than minipig ones. Discussion Single alveolus model is useful to study macrophages-bacilli interaction observing two possible outcomes (three zones of parameters can be distinguished). In the bistability both final states and the probabilities of finishing in each one can be computed. The parameters that reduce bacillary load and exudative scenario probability are identified (#μ, #λT, "ξ, "κ, "ν and "τm) these parameters can be related with biological processes. Then, biological processes that reduce bacillary load and exudative scenario probability are a decrease of bacilli reproduction, an increase in T cells lifespan, an increase in activated macrophages lifespan, a faster activated macrophages killing rate, a faster T cells recruitment speed or a faster macrophages activation rate. Using baseline parameters infection can not be controlled and bacilli resist inflammatory and immune response. The encapsulation process is needed (secondary lobule model) to stop PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 16 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection bacilli proliferation. This process controls lesion growth when immune response is not enough as seen in [7] and in the single alveolus model. The bacilli growth is not stopped but the affected alveoli are controlled not allowing to infect surrounding ones. The modeling of such encapsulation process is closely related with the description of the growth of the granuloma. There have been several attempts to model such growing basically based in agent based models for the cells and their motion among aveoulus of the same sec- ondary lobule [13, 17] and in the other hand reaction diffusion models for the number of cells inside the secondary lobule [43]. Our model also describes the growing of the granuloma by the use of reaction-diffusion equations, however, we incorporate in the model the saturation of such growing due to the interaction with the interlobular connective tissue by the descrip- tion of the dynamics of the fibroblast. As consequence of the approach, lesions size is affected by its position on the lung, in partic- ular, the distance from its focus to the nearest pulmonary membrane. Other two factors that determine its size is the small diffusion velocity of the elements that is reduced by fibroblasts, collagen and other molecules present in alveoli (encapsulation velocity) and the facility that infected macrophages can diffuse to surrounding alveoli (lesion growth). Secondary lobule size and shape is important to control lesions before, therefore to control infection and not evolve to ATB. We have seen that smaller secondary lobule (smaller lungs) present smaller lesions and TB is better controlled. This could be a reason to justify the lower incidence in chil- dren than in adults observed in [44]. Volume and radius evolution of the lesions can be successfully adjusted to a generalized logistic equation as predicted by the bubble model employed in some previous empirical approaches for the growth of the granuloma [17, 18]. Further work Next step is to build a computational lung that can reproduce TB disease that includes: bron- chial tree structure [18] with a secondary lobules structure that includes pulmonary mem- branes and the model described before. The details of the model can be refined including some intermediate steps that were simpli- fied in the processes of macrophages clearance or immune activation. To study post lesion formation scenario it must be considered that when encapsulation is completed foamy macro- phages can not travel outside the lesion, then process 14 (intracellular bacilli drainage) is inhib- ited in that lesion. Another important point that may be taken into account is how ventilation and blood flow affect the different parts of the lung. In particular, secondary pulmonary lobules are in differ- ent conditions depending on the lung part: ventilation, blood flow, shape and size are different. These conditions may change some parameters and end in different lesions outcome. Lesions much bigger than a single secondary lobule are obtained in active tuberculosis. In order to model the break of the pulmonary membrane, some other mechanical processes could be considered. Supporting information S1 Video. Encapsulation of a tuberculosis lesion in a 3D polygon secondary lobule. In green fibroblast that surrounds secondary lobule and encapsulates the lesion. In blue cells that form the lesion (macrophages, neutrophils, T cells. . .). Lesion grows till it contacts with inter- lobular connective tissue and it is encapsulated. (AVI) PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 17 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection Author Contributions Conceptualization: Martı́ Català, Clara Prats, Daniel López, Pere-Joan Cardona, Sergio Alonso. Funding acquisition: Martı́ Català, Clara Prats, Daniel López, Pere-Joan Cardona, Sergio Alonso. Investigation: Martı́ Català, Clara Prats, Daniel López, Pere-Joan Cardona, Sergio Alonso. Methodology: Martı́ Català, Clara Prats, Daniel López, Pere-Joan Cardona, Sergio Alonso. Software: Martı́ Català. Writing – review & editing: Martı́ Català, Clara Prats, Daniel López, Pere-Joan Cardona, Ser- gio Alonso. References 1. World Health Organization. Global Tuberculosis report. World Health Organization; 2018. Available from: https://www.who.int/tb/publications/global_report/gtbr2018_main_text_28Feb2019.pdf?ua = 1. 2. Public health emergency: WHO Global Tuberculosis Programme. TB: a global emergency, WHO report on the TB epidemic. World Health Organization; 1994. Available from: https://apps.who.int/iris/handle/ 10665/58749. 3. World Health Organization. The end TB strategy. World Health Organization; 2015. Available from: https://www.who.int/tb/End_TB_brochure.pdf. 4. Cardona PJ, Català M, Arch M, Arias L, Alonso S, Cardona P, et al. Can systems immunology lead tuberculosis eradication? Current Opinion in Systems Biology. 2018; 12:53—60. https://doi.org/10. 1016/j.coisb.2018.10.004 5. Bermudez L, Danelishvili L, Early J. Mycobacteria and Macrophage Apoptosis: Complex Struggle for Survival. Microbe Magazine. 2006; 1(8). https://doi.org/10.1128/microbe.1.372.1 6. Cardona PJ. Patogénesis de la tuberculosis y otras micobacteriosis. Enfermedades Infecciosas y Microbiologı́a Clı́nica. 2018; 36(1):38–46. https://doi.org/10.1016/j.eimc.2017.10.015 7. Gil O, Dı́az I, Vilaplana C, Tapia G, Dı́az J, Fort M, et al. Granuloma encapsulation is a key factor for containing tuberculosis infection in minipigs. PLoS One. 2010; 5(4):e10030. https://doi.org/10.1371/ journal.pone.0010030 8. Cardona PJ, Prats C. The small breathing amplitude at the upper lobes favors the attraction of polymor- phonuclear neutrophils to Mycobacterium tuberculosis lesions and helps to understand the evolution toward active disease in an individual-based model. Frontiers in microbiology. 2016; 7:354. 9. Cardona PJ. The key role of exudative lesions and their encapsulation: lessons learned from the pathol- ogy of human pulmonary tuberculosis. Frontiers in microbiology. 2015; 6:612. https://doi.org/10.3389/ fmicb.2015.00612 10. Cardona PJ. What we have learned and what we have missed in tuberculosis pathophysiology for a new vaccine design: searching for the “Pink Swan”. Frontiers in immunology. 2017; 8:556. 11. Verrall AJ, G Netea M, Alisjahbana B, Hill PC, van Crevel R. Early clearance of Mycobacterium tubercu- losis: a new frontier in prevention. Immunology. 2014; 141(4):506–513. https://doi.org/10.1111/imm. 12223 12. Haefner JW. Modeling Biological Systems. Springer; 2005. 13. Segovia-Juarez JL, Ganguli S, Kirschner D. Identifying control mechanisms of granuloma formation during M. tuberculosis infection using an agent-based model. Journal of Theoretical Biology. 2004; 231 (3):357—376. https://doi.org/10.1016/j.jtbi.2004.06.031 14. Marino S, Kirschner D. A multi-compartment hybrid computational model predicts key roles for dendritic cells in tuberculosis infection. Computation. 2016; 4(4):39. https://doi.org/10.3390/ computation4040039 15. Warsinske HC, DiFazio RM, Linderman JJ, Flynn JL, Kirschner DE. Identifying mechanisms driving for- mation of granuloma-associated fibrosis during Mycobacterium tuberculosis infection. Journal of theo- retical biology. 2017; 429:1–17. https://doi.org/10.1016/j.jtbi.2017.06.017 16. Ziraldo C, Gong C, Kirschner DE, Linderman JJ. Strategic priming with multiple antigens can yield mem- ory cell phenotypes optimized for Infection with Mycobacterium tuberculosis: a computational study. Frontiers in microbiology. 2016; 6:1477. https://doi.org/10.3389/fmicb.2015.01477 PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 18 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection 17. Prats C, Vilaplana C, Valls J, Marzo E, Cardona PJ, López D. Local Inflammation, Dissemination and Coalescence of Lesions Are Key for the Progression toward Active Tuberculosis: The Bubble Model. Frontiers in Microbiology. 2016; 7:33. https://doi.org/10.3389/fmicb.2016.00033 18. Català M, Bechini J, Tenesa M, Pérez R, Moya M, Vilaplana C, et al. Modelling the dynamics of tubercu- losis lesions in a virtual lung: role of bronchial tree on endogenous reinfection. PLOS Computational Biology. 2020 May 20; 16(5):e1007772. https://doi.org/10.1371/journal.pcbi.1007772 19. Ibargüen-Mondragón E, Esteva L, Burbano-Rosero EM. Mathematical model for the growth of Myco- bacterium tuberculosis in the granuloma. Mathematical Biosciences and Engineering. 2018; 15. http:// dx.doi.org/10.3934/mbe.2018018. 20. Marzo E, Vilaplana C, Tapia G, Diaz J, Garcia V, Cardona PJ. Damaging role of neutrophilic infiltration in a mouse model of progressive tuberculosis. Tuberculosis. 2014; 94(1):55—64. https://doi.org/10. 1016/j.tube.2013.09.004 21. Russell DG, Cardona PJ, Kim MJ, Allain S, Altare F. Foamy macrophages and the progression of the human tuberculosis granuloma. Nature immunology. 2009; 10(9):943–948. https://doi.org/10.1038/ni. 1781 22. Cáceres N, Tapia G, Ojanguren I, Altare F, Gil O, Pinto S, et al. Evolution of foamy macrophages in the pulmonary granulomas of experimental tuberculosis models. Tuberculosis. 2009; 89(2):175–182. https://doi.org/10.1016/j.tube.2008.11.001 23. Cardona PJ. A dynamic reinfection hypothesis of latent tuberculosis infection. Infection. 2009; 37(2):80. https://doi.org/10.1007/s15010-008-8087-y 24. Tak T, Tesselaar K, Pillay J, Borghans JA, Koenderman L. What’s your age again? Determination of human neutrophil half-lives revisited. Journal of leukocyte biology. 2013; 94(4):595–601. https://doi.org/ 10.1189/jlb.1112571 25. Ramos-Kichik V, Mondragón-Flores R, Mondragón-Castelán M, Gonzalez-Pozos S, Muñiz-Hernandez S, Rojas-Espinosa O, et al. Neutrophil extracellular traps are induced by Mycobacterium tuberculosis. Tuberculosis. 2009; 89(1):29–37. https://doi.org/10.1016/j.tube.2008.09.009 26. Bru A, Cardona PJ. Mathematical modeling of tuberculosis bacillary counts and cellular populations in the organs of infected mice. PLoS One. 2010; 5(9):e12985. https://doi.org/10.1371/journal.pone.0012985 27. Mahamed D, Boulle M, Ganga Y, Mc Arthur C, Skroch S, Oom L, et al. Intracellular growth of Mycobac- terium tuberculosis after macrophage cell death leads to serial killing of host cells. Elife. 2017; 6: e22028. https://doi.org/10.7554/eLife.28205 28. Riley LK, Rupert J. Evaluation of patients with leukocytosis. American family physician. 2015; 92 (11):1004–1011. 29. Kramer A, Schwebke I, Kampf G. How long do nosocomial pathogens persist on inanimate surfaces? A systematic review. BMC infectious diseases. 2006; 6(1):130. https://doi.org/10.1186/1471-2334-6-130 30. Murphy J, Summer R, Wilson AA, Kotton DN, Fine A. The Prolonged Life-Span of Alveolar Macro- phages. American Journal of Respiratory Cell and Molecular Biology. 2008; 38(4):380–385. https://doi. org/10.1165/rcmb.2007-0224RC 31. Gammack D, Doering C, Kirschner D. Macrophage response to Mycobacterium tuberculosis infection. Journal of mathematical biology. 2004; 48:218–42. https://doi.org/10.1007/s00285-003-0232-8 32. Krombach F, Münzing S, Allmeling AM, Gerlach JT, Behr J, Dörger M. Cell size of alveolar macro- phages: an interspecies comparison. Environmental Health Perspectives. 1997; 105(suppl 5):1261– 1263. https://doi.org/10.2307/3433544 33. Ting-Beall H, Needham D, Hochmuth R. Volume and osmotic properties of human neutrophils. Blood. 1993; 81(10):2774–2780. https://doi.org/10.1182/blood.V81.10.2774.2774 34. Chapman EH, Kurec AS, Davey FR. Cell volumes of normal and malignant mononuclear cells. Journal of Clinical Pathology. 1981; 34(10):1083–1090. 35. Saladin K. Human anatomy. 3rd ed. McGraw-Hill; 2011. 36. Todar K. Online textbook of bacteriology. University of Wisconsin-Madison Department of Bacteriology. 2008;. 37. Fürth R. Die brownsche bewegung bei berücksichtigung einer persistenz der bewegungsrichtung. mit anwendungen auf die bewegung lebender infusorien. Zeitschrift für Physik A Hadrons and Nuclei. 1920; 2(3):244–256. 38. Wheeler AP, Wells CM, Smith SD, Vega FM, Henderson RB, Tybulewicz VL, et al. Rac1 and Rac2 reg- ulate macrophage morphology but are not essential for migration. Journal of Cell Science. 2006; 119 (13):2749–2757. https://doi.org/10.1242/jcs.03024 39. Gillespie DT. Approximate accelerated stochastic simulation of chemically reacting systems. The Jour- nal of chemical physics. 2001; 115(4):1716–1733. https://doi.org/10.1063/1.1378322 PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 19 / 20 PLOS ONE Reaction-diffusion model to understand granulomas formation in secondary lobule during tuberculosis infection 40. Marino S, Hogue IB, Ray CJ, Kirschner DE. A methodology for performing global uncertainty and sensi- tivity analysis in systems biology. Journal of theoretical biology. 2008; 254(1):178–196. https://doi.org/ 10.1016/j.jtbi.2008.04.011 41. McKay MD, Beckman RJ, Conover WJ. Comparison of three methods for selecting values of input vari- ables in the analysis of output from a computer code. Technometrics. 1979; 21(2):239–245. https://doi. org/10.1080/00401706.1979.10489755 42. Richards F. A flexible growth function for empirical use. Journal of experimental Botany. 1959; 10 (2):290–301. https://doi.org/10.1093/jxb/10.2.290 43. Hao W, Schlesinger LS, Friedman A. Modeling granulomas in response to infection in the lung. PLoS One. 2016 Mar 17; 11(3):e0148738. https://doi.org/10.1371/journal.pone.0148738 44. Seddon JA, Shingadia D. Epidemiology and disease burden of tuberculosis in children: a global per- spective. Infection and Drug Resistance. 2014; 7:153. PLOS ONE | https://doi.org/10.1371/journal.pone.0239289 September 16, 2020 20 / 20