Modelling the dynamics of tuberculosis lesions in a virtual lung: role of the bronchial tree in endogenous reinfection

RESEARCH ARTICLE Modelling the dynamics of tuberculosis lesions in a virtual lung: Role of the bronchial tree in endogenous reinfection Martı ´Català 1,2 , Jordi Bechini 3 , Montserrat Tenesa 3 , Ricardo Pe ´rez 3 , Mariano Moya 3 , Cristina VilaplanaID 4,5 , Joaquim VallsID 2 , Sergio AlonsoID 2 , Daniel Lo ´pez 2 , PereJoan CardonaID 1,4,5 , Clara PratsID 2 * 1Comparative Medicine and Bioimage Centre of Catalonia (CMCiB), Fundacio ´Institut d’Investigacio ´en Ciències de la Salut Germans Trias i Pujol, Badalona, Catalonia, Spain, 2Departament de Fı ´sica, Universitat Politècnica de Catalunya, Castelldefels, Barcelona, Catalonia, Spain, 3Servei de Radiodiagnòstic, Hospital Universitari Germans Trias i Pujol, Badalona, Catalonia, Spain, 4Experimental Tuberculosis Unit, Fundacio ´ Institut d’Investigacio ´en Ciències de la Salut Germans Trias i Pujol, Universitat Autònoma de Barcelona, Can Ruti Campus, Edifici Mar, Badalona, Catalonia, Spain, 5Centro de Investigacio ´n Biome ´dica en Red de Enfermedades Respiratorias, Madrid, Spain *[email protected] Abstract Tuberculosis (TB) is an infectious disease that still causes more than 1.5 million deaths annually. The World Health Organization estimates that around 30% of the world’s population is latently infected. However, the mechanisms responsible for 10% of this reserve (i.e., of the latently infected population) developing an active disease are not fully understood, yet. The dynamic hypothesis suggests that endogenous reinfection has an important role in maintaining latent infection. In order to examine this hypothesis for falsifiability, an agentbased model of growth, merging, and proliferation of TB lesions was implemented in a computational bronchial tree, built with an iterative algorithm for the generation of bronchial bifurcations and tubes applied inside a virtual 3D pulmonary surface. The computational model was fed and parameterized with computed tomography (CT) experimental data from 5 latently infected minipigs. First, we used CT images to reconstruct the virtual pulmonary surfaces where bronchial trees are built. Then, CT data about TB lesion’ size and location to each minipig were used in the parameterization process. The model’s outcome provides spatial and size distributions of TB lesions that successfully reproduced experimental data, thus reinforcing the role of the bronchial tree as the spatial structure triggering endogenous reinfection. A sensitivity analysis of the model shows that the final number of lesions is strongly related with the endogenous reinfection frequency and maximum growth rate of the lesions, while their mean diameter mainly depends on the spatial spreading of new lesions and the maximum radius. Finally, the model was used as an in silico experimental platform to explore the transition from latent infection to active disease, identifying two main triggering factors: a high inflammatory response and the combination of a moderate inflammatory response with a small breathing amplitude. PLOS COMPUTATIONAL BIOLOGY PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 1 / 25 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: CatalàM, Bechini J, Tenesa M, Pe ´rez R, Moya M, Vilaplana C, et al. (2020) Modelling the dynamics of tuberculosis lesions in a virtual lung: Role of the bronchial tree in endogenous reinfection. PLoS Comput Biol 16(5): e1007772. https://doi.org/10.1371/journal.pcbi.1007772 Editor: Dominik Wodarz, University of California Irvine, UNITED STATES Received: June 5, 2019 Accepted: March 4, 2020 Published: May 20, 2020 Copyright: ©2020 Catalàet al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: All relevant data are within the manuscript and its Supporting Information files. Funding: CP, PJC and MC received funding from “la Caixa” Foundation (ID 100010434), under agreement LCF/PR/GN17/50300003; CV, PJC, MM, JB, MT and RP received funding from Agència de Gestio ´d’Ajuts Universitaris i de Recerca (AGAUR), Grup Unitat de Tuberculosi Experimental, 2017SGR-500; CP, DL, SA, MC, JV received funding from Ministerio de Ciencia, Innovacio ´n y Author summary Tuberculosis is, even today, among the 10 main causes of death in the world. Despite the effectiveness of current strategies to fight the disease and those that are under development, the huge reservoir of latently infected individuals is a big hindrance in its eradication. One of the challenges inherent in this problem is that the mechanisms that cause latent infection to evolve towards active disease are not fully understood. Why will 90% of infected individuals never develop an active disease? In other words, what are the main factors that trigger an active disease in 10% of cases? We have focused our efforts on understanding the mechanisms that allow keeping infection latent, especially those related with endogenous reinfection. Since it is supposed to occur through the bronchial tree, we have designed a 3D computational model that mimics this structure, in which we have implemented an agent-based model of lesion growth and proliferation. Our results were contrasted with computed tomography measurements in latently infected minipigs, providing successful results that reinforce the essential role of endogenous reinfection through the bronchial tree in keeping infection latent. Introduction Tuberculosis (TB) is an infectious disease that in 2017 killed more than 1.6 million people. Mycobacterium tuberculosis (Mtb) causes TB, and this bacterium is the individual agent causing the highest mortality worldwide [1]. The World Health Organization (WHO) estimates that 25 to 30% of the population worldwide is infected with Mtb, and that around 10% of infected people will develop active tuberculosis (ATB) in a few years’ time [2], although these percentages are being questioned and re-visited by recent studies [3]. WHO also estimates that 10 million humans developed ATB in 2017 [2]. TB infection starts at a pulmonary alveolus when Mtb is phagocyted by an alveolar macrophage (AM). The Mtb resists bactericidal mechanisms induced by AM and replicates inside the macrophage [4]. Under proper in vitro conditions Mtb replicates once a day [5]. When the intracellular bacterial load overcomes the AM’s maximum tolerability, macrophage necrosis is triggered, thereby returning bacilli to the extracellular milieu. These bacilli are phagocyted by other AMs and the cycle begins again giving rise to a further increase in bacilli. The further inclusion of more AMs fails to control bacillary growth. The death of AMs triggers a local inflammatory response first, and then a specific immune response, which finally controls the infection. The end of the progressive infection leaves an encapsulated TB lesion [6]. According to the dynamic hypothesis of Cardona [7], there is a certain probability that a few bacilli will escape from the lesion, mostly inside a foamy macrophage, and start a new infection in another alveolus. This process is assumed to occur through the bronchial tree, and it is what we denote as endogenous reinfection. This includes not only new infections generated from the initial infection site, but also those originated in successive infection foci (Fig 1). An Mtb infection may be completely cleared by the organism [8], it may enter into a latent state in which the host is infected but not sick and cannot infect other people, called Latent Tuberculosis Infection (LTBI), or, if the immune and inflammatory responses are not well balanced, the host may develop Active Tuberculosis disease (ATB). Systems biology and computational models are fruitful tools for increasing understanding of the processes involved in TB [1]. Recently, different models have been useful in identifying several TB key factors [9–13]. In particular, the Bubble model suggests that the coalescence of closed lesions is the main mechanism for the growth of lesions in animals that progress to PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 2 / 25 Universidades and FEDER, with the project PGC2018-095456-B-I00; CVM has a Miguel Servet II contract funded by the Instituto Carlos III (ISCIII, CPII18/00031). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist. ATB [14]. This model successfully explains experimental observations in mice [15]. The Bubble model assumes a generalised logistic (Richard’s curve) [16] growth of lesions, driven by the inflammatory and immune responses, with their proliferation according to the endogenous reinfection theory, and a merging between neighbouring lesions when they are close enough. The model successfully reproduced ATB observed in C3HeB/FeJ mice, demonstrating the importance of local inflammation, lesion proliferation, and coalescence in the triggering of active disease. These results are relevant for mouse models; however, they are not easily extrapolated to humans, because of the differences between the structure of the lungs in the two species, in addition to the well-known differences in immune systems and encapsulation capacity. Actually, the structure of the lungs may play an important role in the infection dynamics of TB. On the one hand, endogenous reinfection occurs mainly through the bronchial tree, and mice have much simpler pulmonary structure than humans, as no secondary lobular structure is found in mice (they have little or none interlobular septae) [17]. On the other hand, the encapsulation of lesions is driven by fibroblasts and fibrin from pulmonary membranes like intralobular septae. Nevertheless, mice do not possess intralobular septae and lesion encapsulation is not possible except close to the pleura [18]. In addition, the immune response of C3HeB/FeJ mice is much less effective than that of humans. In humans, a balanced Th1 immune response usually takes place after TB infection [19]. As mice’s immune response is not strong and encapsulation is not possible, these animal models cannot develop an LTBI situation and all experimental observations show ATB cases [15]. Although pigs and humans share a great deal of anatomy and physiology, researchers rarely employ pigs as in vivo models for TB. Yet their immune system and lung structure are particularly close to the corresponding system and structure in humans. Thus, TB development in Fig 1. Main features of dynamic hypothesis. Schematic representation of initial infection (A, white arrow) lesions’ growth (A– F), endogenous reinfection (C–E, yellow arrows) and coalescence of neighbouring lesions (F). https://doi.org/10.1371/journal.pcbi.1007772.g001 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 3 / 25 pigs is more similar to that in humans than in mouse models [19,20]. Minipigs are a genetically selected species, which is more convenient than other pigs for experiments in a lab, mainly for size reasons. Experimental results in TB in minipigs resemble pathological findings described in human [21–23]. In this study we aim to adapt and implement the Bubble model in a virtual bronchial tree in order to understand the maintenance of LTBI in minipigs. In particular, we want to test the falsifiability of the dynamic hypothesis of Cardona [7] that explains this maintenance, as well as to obtain some orders of magnitude of its dynamics. We use experimental minipig TB data to tune the model [23]. With the new model we perform several in silico experiments, which successfully reproduce experimental observations, and, furthermore, permit us to systematically explore the transition between ATB and LTBI. In Materials and Methods we describe the CT experimental data, as well as the two submodels used in simulations, which correspond to the computational lung and the revised Bubble model. We finish this section providing details of the model’s implementation and the methodology used for its parameterization and sensitivity analysis. Results’ section starts with an analysis of the computational lung obtained. Then, it provides the results of the model’s fitting to experimental data and the sensitivity analysis. Subsequent simulation series are used to explore the effect of the initial configuration and to test the transition between LTBI and ATB in minipigs. Finally, the conclusions for this study and their implications in testing the falsifiability of the dynamic hypothesis are drawn. Materials and methods Computer tomography measurements of LTBI in minipigs Six female specific pathogen-free (spf) minipigs were intratracheally infected by H37Rv Pasteur strains of Mtb (10 3 CFU) under sedation [23]. They were euthanized twelve weeks postinfection, without having received any TB treatment. None of them had developed ATB symptoms. Broncho-pulmonary pieces were obtained from all minipigs and analysed with multidetector computed tomography scan (CT), a GE LightSpeed VCT with high image resolution (64-slice). The morphoanatomical study was carried out with volume rendering software. For each of the minipigs, we recorded the number of lesions, their location and size, Hounsfield units, and distance to closest pleura. One of the analysed minipigs did not present any lesions; it was considered non-infected for technical reasons and excluded from the subsequent analysis. The 5 infected animals showed 165 lesions in total, 33 ±22 per minipig. These lesions are shown in Fig 2. The mean diameter of the lesions was 1.3 ±0.2 mm. Lesions were located in each minipig with a mean dispersion of 16 ±4 mm. The number of lesions and their positions were used to train the computational model [24]. Ethics statement. All ethical requirements were followed according to Directive 201/63/ EU, and the protocol and procedures of the study were approved by the corresponding ethical committee on animal welfare and the Catalan Government (Permit number: 5796). All animals were euthanized at week 12 post-infection by intravenous injection of sodium pentobarbital. Computational bronchial tree The main novelty of this modelling approach is the use of an explicit 3D space that resembles a pulmonary bronchial tree. The design of this explicit space requires the building of a computational bronchial tree inside a certain pulmonary volume, limited by the external surface. The PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 4 / 25 geometrical information necessary for this model can be obtained from pulmonary CT images of the studied minipigs. Accordingly, the computational bronchial tree consists of two models: (1) an empirical model for the external lung’ surface that limits pulmonary volume, and (2) an artificial iterative model of bifurcations to build a bronchial tree inside this surface. This model is deterministic, since surfaces are obtained from experimental CT measurements and the iterative model does not incorporate randomness. Empirical model for the lung’s external surface. An empirical surface model is built using CT scan data from one of the minipigs, randomly chosen. This representative surface is subsequently re-scaled according to the dimensions of others minipigs’ lungs, giving rise to 5 computational surfaces that can be used to buildi the 5 different bronchial trees. In order to build the representative pulmonary surface, we use three images of the three planes, i.e., coronal plane, sagittal plane, and axial plane. From these images, the contour line is extracted, keeping the carina (i.e., the bifurcation point of the trachea where it divides into the two main bronchi) position for purposes of reconstruction (Fig 3). The 3D reconstruction from contour lines is carried out with Matlab. All contours are normalized to 1 in order to be subsequently re-scaled with the specific dimensions of each of the 5 minipigs’ lungs on the reconstruction process. Finally, the left lung is slightly rotated (5˚) so that the inter-pulmonary space is reduced. All the lesions observed experimentally are located in the coordinate system defined by the reconstructions, taking the carina as the reference point, in order to check if they are located inside the obtained computational pulmonary surfaces. As a result, 95% of the detected lesions are inside the computed surface or in contact with pleura. Iterative model for the bronchial tree. A bronchial tree of the conductive zone is built inside the computational pulmonary surface with an algorithm based on previous work on the human bronchial tree [25,26]. Starting at the trachea, a set of iterative rules govern the successive bifurcations. The bronchial tree of the minipigs is assumed to be morphologically equivalent to the human one, but with smaller dimensions [27]. Fig 2. Summary of experimental results. Left: CT image reconstruction of a minipig’s pulmonary surface. Right: 3D representation of location and size of all minipig lesions; each colour is for a different minipig (red, magenta, blue, yellow, and green). https://doi.org/10.1371/journal.pcbi.1007772.g002 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 5 / 25 Our algorithm assumes that all the divisions are bifurcations, i.e., they occur in a dichotomous way. The resulting three branches involved in a bifurcation are coplanar, and the plane that contains each bifurcation is called a bifurcation plane. The divisions are assumed to occur in successive perpendicular planes, i.e., right-left, anterior-posterior, and upper-lower. Therefore, each bifurcation plane is perpendicular to the previous plane. The first division starts at the carina and directs the new branches into the right and left lungs. When a certain conducting airway 0 divides into conducting airways 1 and 2, the flow conservation (Q 0 = Q 1 + Q 2 ) together with Murray’s law (Q = C �d 3 ) [28] leads to the following relation between their diameters, d i : d0 3¼d1 3þd2 3ð1Þ Florens et al. [29] derived a ratio of 3 for the length of a branch (l i ) and its diameter (d i ) for most of the bronchial trees: li¼3dið2Þ We also studied this relation using Rozanek and Roubik’s experimental data [30], obtaining a proportionality constant of 3.07 and a goodness of fit of R 2 = 0.98. This analysis is shown in the Supplementary material section 3 (Fig B in S1 File). Fig 3. Normalized contour lines obtained from CT-scan images. Contour lines were obtained from CT-scan images and used for a 3D computational reconstruction of the pulmonary surface. Trachea division point (carina) is marked for purposes of reconstruction. (A) Sagittal plane CT image. (B) Coronal plane CT image. (C) Axial plane CT image. (D) Sagittal plane outline reconstruction. (E) Coronal plane outline reconstruction. (F) Axial plane outline reconstruction. https://doi.org/10.1371/journal.pcbi.1007772.g003 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 6 / 25 The diameters and angles of each bifurcation depend on the relative volume that each new branch supplies. We define the q i factor as: qi¼Vi V0 ;i¼1;2ð3Þ where V 1 and V 2 are, respectively, the sub-volumes irrigated by conducting airways 1 and 2 after bifurcation, and V 0 is the volume supplied by the branch 0 (before bifurcation). Taking into account Murray’s law, the diameters after the bifurcation are: di¼d0�qi 1=3;i¼1;2ð4Þ Minimizing the work per unit time (associated with friction and to maintain the structure) in bifurcations, the following relations between the angles of the bifurcation and the factor q i are obtained [31]: cos�i¼1 2qi2=31þqi 4=3 ð1qiÞ4=3 � �;i¼1;2ð5Þ The calculation of the ratio q i (Eq 3) cannot be analytically evaluated. Therefore, a grid of equispaced points is created so that the number of points inside each considered volume, N i (i = 0,1,2), is assessed and the ratio is evaluated as: qi¼Ni N0 ;i¼1;2ð6Þ The distance between points is initially fixed at 1 mm [25] and then reduced to 0.2 mm to increase precision and improve results. In Fig 4 there can be seen a diagram of a bifurcation example for q 1 = 0.6. Using Eq 4 it can be determined that the diameter of the daughter branches are d 1 = 0.84�d 0 and d 2 = 0.74�d 0 , respectively. Length is 3 times the diameter of each branch, then: l 0 = 3�d 0 , l 1 = 3�d 1 = 2.53�d 0 = 0.84�l 0 and l 1 = 3�d 1 = 2.21�d 0 = 0.74�l 0 . Bifurcation angles can be computed using Eq 5 as: ϕ 1 = 32˚ and ϕ 2 = 43˚. The initial tests with these equations resulted in a bronchial tree that does not fully occupy the upper part of the pulmonary surface. Therefore, a correction is added to the angle calculation [25]. It consists of an evaluation of the mass centre of the volumes to be supplied, !xMC;Vi, which is projected on the bifurcation plane. We will denote with ψ i the angle between the original branch and the projected mass centre. Experimentally, bifurcations with an angle greater than π/2 were not observed [32]. Therefore, the corrected angle φ i is: φi¼min �iþci 2;p 2 � �;i¼1;2ð7Þ Branches with a diameter lower than 0.5 mm are considered terminals as areas those branches that escape from the pulmonary surface. The model for the building of successive bifurcations is summarized in Table 1. The Bubble model: an update The Bubble model is an agent-based model in which the agents are the lesions. The dynamics of the lesions are driven by three processes: a generalised logistic growth, the reinfection process that permits the generation of new infection focuses, and the coalescence between neighbouring lesions. PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 7 / 25 The Bubble model was originally designed and calibrated to describe the dynamics of tuberculous lesions in mice with an active disease [14]. The spatial structure was not relevant for mice, due to the dimensions and the relatively simple structure of mouse lungs, and taking into consideration the size of lesions for the active disease. This model is updated and implemented as follows. Lesion growth. We model a lesion as a sphere whose spatial position (3D coordinates, in mm), radius (in mm), and age (in days) are variables. Lesions are firstly detected when their radius is r min = 0.075 mm (smaller lesions cannot be identified). This occurs after approximately t min = 14 days from the initial infection. Then when a lesion is created it remains “silent” for 14 days before it is initialized with a r min radius. The model employs a generalised logistic growth of the radius of the lesions as follows: driðtÞ dt ¼vi�ritð Þ� 1riðtÞ rmax � �2 " # ð8Þ where r i is the radius of the lesion, ν i is the parameter that sets the maximum growth rate, and r max is the maximum radius. The parameter ν i is modelled as a Gaussian variable with mean value νand standard deviation ν/3. Therefore, each lesion grows at a slightly different velocity at each time step. From experimental data it is known that around the 28 th day a 2 mm lesion reaches its limit [22]; therefore, νis estimated as ν= 0.3 day -1 (σ ν = 0.1 day -1 ). Fig 4. Bifurcation diagram. Bifurcation of 0 branch into two (1, 2) daughter branches. The cabal ratio for branch 1 is: q 1 = 0.6. Length is 3 times the diameter of each branch as may be seen in Eq 2. Diameter relations are obtained from Eq 4, as d 1 = 0.84�d 0 and d 2 = 0.74�d 0 . Angular values are computed using Eq 5, as ϕ 1 = 32˚ and ϕ 2 = 43˚. https://doi.org/10.1371/journal.pcbi.1007772.g004 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 8 / 25 Lesion proliferation. The multiplication of the number of lesions is caused by endogenous reinfection. In this way, a mother lesion generates new daughter lesions from day 14 to day 28. The original reinfection probability function [14] includes two terms: (1) a linearly increasing term with the radius of the mother, and (2) a linearly decreasing probability with mother lesion age [22]. The second term is slightly modified in order to allow the generation of new daughter lesions from mother lesions older than 28 days, with a small non-zero probability: P tð Þdt ¼r�riðtÞ rmin �eaðai14Þndt;ai�14days ð9Þ where ρ,α, and nare parameters that define the probability profile, a i is the age of the lesion, in days, and r min is the minimum radius at which lesions are identified. Fig 5 shows original [14] and modified (Eq 9) models with α= 0.035 day -n and n= 1.63 with r max = 1 mm and ρ= 0.10 day -1 . The values of αand nare fixed to ensure that the area under the two curves are equal and to minimize the difference between the experimental results and the linear model. The reinfection probability initially grows exponentially with the lesion radius; however, for longer times the curve decays exponentially because of encapsulation and calcification. In the original model for mice, the location of new lesions is selected from a probability which decreases with the distance. In the current version for minipigs the distance is modelled in the same way, but considering the bronchial distance between two terminals instead of the geometric distance between two points. The model does not explicitly assume that the dissemination occurs exclusively through the aerial part, but it could include a possible recirculation through the adjacent circulatory or lymphatic systems as well. Then, possible locations are determined by the bronchial tree terminal positions: P i !jð Þ ¼ ebdij Pjebdij ð10Þ where P(i!j) is the probability that a lesion appears at a terminal jdue to a mother lesion at terminal i, and βis the dispersion parameter that determines the spreading. Fig 6A shows the mean distance of the appearance of new lesions as a function of the dispersion parameter. Fig Table 1. Summary of the model for building the computational bronchial tree of each minipig. The normalized pulmonary surface is specifically re-scaled for each minipig, taking into account the measured dimensions. A bronchial tree is built inside each computational pulmonary surface. The bronchial tree starts at the end of the trachea, taking the trachea diameter and carina location of each minipig as reference. The bronchial tree is a tubular structure, and non-terminal branches split in a dichotomous way. The three branches implied in a division are coplanar (bifurcation plane). In each division, the pulmonary territory is divided in two subregions by the plane that is perpendicular to the bifurcation plane, following the direction of the mother branch. The bifurcation plane of the first division is the vertical one that separates the right and the left lungs. The bifurcation planes are perpendicular from one generation to the following. Once the bifurcation plane is defined, the ratio q i =V i /V 0 is numerically evaluated by a grid of equiespaced points with a precision of 0.2 mm, so that q i =N i /N 0 (N i and N 0 are the number of points contained in each subvolume). The parameters of each bifurcation are evaluated using Eqs 2,4,5and 7. Branches that escape the pulmonary surface and those with a diameter of less than 0.5 mm are considered terminals. https://doi.org/10.1371/journal.pcbi.1007772.t001 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 9 / 25 The lack of reliable information about the initial infection entails the need for a blind assumption. Following the law of parsimony, we assume the simplest initial configuration, and alternative possibilities will be explored later on. We use as initial infection for the simulations a single lesion located at the nearest terminal to the mass centre of all the lesions, given the final distribution shown by CT images. The set of parameters is adjusted to minimize the errors as explained in Materials and Methods. followed by the performance of 500 independent simulations of each minipig’s configuration (i.e., a total of 2500 simulations). The set of parameters obtained after the minimization of the objective functions (Eqs 17–19) is ρ= 0.13 day -1 ,β= 0.08 mm -1 , and r max = 0.68 mm, which corresponds to NLE = 0.089, DE = 0.283, and SE = 0.299. In Fig 9 we show the comparison between the experimental and the simulated outcome distributions of spatial coordinates X, Y, and Z, and of lesion diameters. A two-sample Kolmogorov-Smirnov test of the four pairs of distributions show that there are not significative differences in any of the cases with a signification level of 0.05. Therefore, the distributions resulting from the numerical simulation successfully reproduce the experimental observations. Simulations show that coalescence of lesions is nearly non-existent, on average less than one coalescence per minipig. This result is in agreement with Prats et al. [9] and Marzo et al. [15], who presented coalescence as a mechanism essential to the evolution towards an active Fig 9. Comparison of experimental and computationally obtained distributions of lesion location and size. Computational distributions were obtained considering as initial infection one lesion in the mass centre of the observed lesions. The set of parameters used is: ρ= 0.13 day -1 ,β= 0.08 mm -1 , and r max = 0.68 mm. (A) Coordinate X (Left— Right) histogram comparison. (B) Coordinate Y (Anterior–Posterior) histogram comparison. (C) Coordinate Z (Vertical) histogram comparison. (D) Diameter distribution histogram comparison. https://doi.org/10.1371/journal.pcbi.1007772.g009 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 16 / 25 disease from a latent infection. A lack of coalescences, therefore, would be a control indicator of the latent infection. Furthermore, simulation results show that the process of endogenous reinfection is crucial to understanding how lesions appear at different locations in the lungs. The use of a computational bronchial tree for driving such reinfection produces spatial distributions which resemble the experimental cases. This process has been shown to be a key factor in maintaining a latent infection inside big mammals like minipigs. Sensitivity analysis Table 5 shows the results of the sensitivity analysis, with the minimum value for the ANOVA test between the increased and the decreased parameters. This analysis reveals that the number of lesions is strongly related with ρ,ν, and T max . It was not obvious that the number of lesions would be related with the growth velocity; however, when lesions grow faster there is an increase in the likelihood that the process of endogenous reinfection will generate new lesions. The mean diameter varies with parameters β,r max , and T max . The inflammatory response is the cause of lesion growth, so relations with r max and T max are expected. The results also show that the dispersion parameter, β, slightly affects the mean diameter. A smaller dispersion parameter causes lesions to be closer, thereby increasing the chance of a coalescence event. In fact, as seen in extended sensitivity analysis for the explored parameter space, r max and βare the two parameters that affect the mean diameter value most. An increase in one of these parameters increases mean diameter value. According to this analysis, T max is the only parameter that is not related with the resulting number of coalescences. All other parameters affect the coalescence processes; nevertheless, a counter-intuitive result is that the dispersion parameter is not the most strongly related. One may expect the dispersion parameter to be the parameter that would affect the coalescence process most because it is the one that determines coalescence spreading, and then determines the distance at which new lesions appear. The sensitivity analysis evidently depends on the initial set of parameters. Nevertheless, we have employed other sets of parameters to carry out sensitivity analyses, providing equivalent results. These analyses are shown in the Supplementary material section 1 (Table A and Table B). Analysing the effect of initial conditions After confirming that the model with the simplest assumption for the initial conditions is good enough to explain the experimental results, we explored the possibility of improving the Table 5. Sensitivity analysis for the set of parameters: S = {ρ,β,r max } = {0.12 day -1 , 0.08 mm -1 , 0.68 mm}. Input parameters βr max ρνT max Outcome variables Number of lesions 0.137 0.092 <0.001�� <0.001�� <0.001�� Mean diameter 0.003�� <0.001�� 0.475 0.106 <0.001�� Dispersion <0.001�� 0.613 0.089 <0.001�� <0.001�� Coalescences 0.024�0.004�� <0.001�� <0.001�� 0.754 Each value is the minimum of the p-values from ANOVA test comparing 500 runs of the original set of parameters and the set of parameters where one parameter is increased or decreased by 10%. Numbers marked with �present statistically significant differences with p<0.05 and numbers marked with �� with p<0.01. https://doi.org/10.1371/journal.pcbi.1007772.t005 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 17 / 25 agreement between the model and experimental measurements by testing different initial distributions of lesions. Table 6 shows the 12 initial configurations analysed, in addition to the previous one. We explore the choice of one or more lesions from CT data as the initial infection using location, size, and density criteria, as well as different random choices. For each initial distribution, the set of parameters is adjusted to minimize the objective functions as explained in Materials and Methods. Table 6 shows the parameter values that minimize errors as well as the corresponding values for each of the explored initial distributions. The results of this analysis, shown in Table 6, do not provide a conclusive criterion for distinguishing those lesions that belonged to the initial infection. Nevertheless, they corroborate that the final lesion distribution is strongly related with the initial infection distribution, since the objective function that is most affected is that of spatial error (SE). The distributions that assume as initial infection one or two random experimental lesions provide better spatial error objective function results; however, they give rise to larger values for DE. In some cases, the value of the spatial parameter that minimizes the spatial error function is β= 0 mm -1 . This value permits a macrophage to travel to all other terminals not taking into account the distance through the bronchial tree. In consequence, the final distribution becomes one that follows all terminal spatial distributions and is not related with the particular initial distribution. Therefore, in these cases the initial distribution is related neither with density nor diameter. From latent to active tuberculosis: In silico experiments Mathematically we define a case of active disease as one with numerical simulations providing a lesion larger than 1 cm in diameter [35]. The model is designed to reproduce experimental results from latent tuberculosis in minipigs. Therefore, no trigger of disease is observed in any of simulations with the fitted parameters. The following set of simulations is designed to explore the parameter space, looking for those zones leading to active disease. Table 6. Initial configurations explored with the model using experimental data. Initial infection Parameters set Errors ρ (day -1 ) β (mm -1 ) r max (mm) NLE DE SE SE/SE control Mass centre (control) 0.134 0.071 0.67 0.015 0.25 0.30 1.00 Coordinate centre 0.134 0.070 0.66 0.026 0.23 0.34 1.14 Biggest lesion 0.129 0 0.68 0.006 0.18 0.74 2.45 Two biggest lesions 0.102 0 0.68 0.015 0.18 0.73 2.42 30% biggest lesions 0.084 0 0.69 0.009 0.17 0.77 2.58 Densest lesion 0.123 0 0.69 0.003 0.19 0.73 2.43 Two densest lesions 0.110 0 0.68 0.026 0.17 0.73 2.44 30% densest lesions 0.084 0 0.68 0.004 0.16 0.78 2.58 Density>150HU 0.078 0 0.68 0.020 0.17 0.79 2.62 One random lesion 0.227 0.170 0.54 0.027 0.63 0.25 0.82 Two random lesions 0.122 0.129 0.58 0.008 0.35 0.23 0.77 One random terminal 0.212 0.160 0.55 0.011 0.59 0.71 2.36 Two random terminals 0.129 0.150 0.59 0.016 0.40 0.71 2.37 The first column shows the criteria for choosing which of the measured lesions were assumed as initial infection. The following columns show the parameter values that minimized errors (ρ,β,r max ) and the values of the three errors obtained (NLE,DE, and SE). The last column compares the spatial error objective function (SE) with the one from the control simulation. HU: Houndsfield units. https://doi.org/10.1371/journal.pcbi.1007772.t006 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 18 / 25 The parameter space is delimited by β2[0, 0.2] mm -1 ,r max 2[1, 5] mm, and ρ2[0.02, 0.2] day -1 . We used equidistant points, 11 for β, 10 for r max , and 4 for ρ. We explored a total of 440 points. We ran 2500 simulations for each point of this parameters space, and 500 for each minipig virtual lung. The initial infection configuration is set as the control (i.e., one lesion in the mass centre of the measured lesions’ distribution). Finally, we define an Active Disease Index as the frequency of active cases among the total number of in silico experiments for each point of the parameter space. We show in Fig 10 the obtained Active Disease Index for each of the explored points. This index increases with any of the 3 parameters, r max ,β, and ρ. Three different zones of the parameters space can be defined: a latent zone, where most of the simulations result in a latent infection (green colour in Fig 10), an active disease zone, where most of the simulations derive into an active disease (red colour in Fig 10), and a transition zone, where both dynamics are possible (from light green to dark orange in Fig 10). The parameter r max is related with the effective inflammatory response in a broad sense; the greater the effect of the inflammatory response, the bigger the lesions and the greater the likelihood of developing an active disease. This result is in agreement with experimental observations [15] and with other biology systems approaches [14]. Of course, the effective dynamics of inflammatory response can be modulated by local properties such as oxygen concentration or macrophages’ availability, among others, which are not explicitly considered by the model. Parameter βis the dispersion parameter; a higher value of βcorresponds to lower dispersion inside the lung, which can be a consequence of a lower breathing amplitude. Therefore, a Fig 10. Active Disease Index for different sets of parameters. Exploration of parameters space (r max ,βand ρ) to see the fraction of in silico experiments that present an active TB disease. The colour is proportional to this frequency; green colour means most of the cases remained latent, red colour means that most of the cases derived into an active disease, and intermediate colours mean that both dynamics are possible. https://doi.org/10.1371/journal.pcbi.1007772.g010 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 19 / 25 low breathing amplitude appears again as a possible cause for the appearance of big lesions, as was previously described in the literature [9]. It has to be taken into account that breathing amplitude changes from one lobe to another one, i.e., it is wider in lower lobes. Therefore, the conditions that facilitate the development of an ATB would vary from one lobe to another. Our model shows that active disease can be triggered by a high inflammatory response or due to a moderate inflammatory response combined with a small breathing amplitude. Nevertheless, at this point the virtual lungs are considered to be homogeneous, i.e., we are not considering variations of parameters along the lungs’ structure. After this exploration, the sensitivity analysis is extended to check if sensitivity of the model depends on the set of parameters used. In addition to the default parameters that belong to the latent infection zone, we chose two combinations of parameters: one representative of the active disease zone and the other of the transition zone; see Fig 10. A latent TB parameter set shows that the space parameter, β, determines mostly the lesions dispersion; however, a parameter set of the transition zone (r max = 6.5 mm, β= 0.1 mm -1 ,ρ= 0.05 day -1 ) shows that the space parameter affects mostly the number of lesions, mean diameter, and coalescences, and does not affect the dispersion of the lesions. This can be seen in Supplementary material, section 1. Discussion Limitations and further work In this modelling approach we have followed the law of parsimony (Occam’s razor) [36], trying to find a simple solution for complex problems such as TB infection dynamics in lungs. The level of complexity was chosen according to the questions to be addressed. This method was developed specifically to submit the main assumptions of the dynamic hypothesis to falsifiability testing. Therefore, the current model includes the most important steps of TB infection evolution suggested by this hypothesis: endogenous reinfection, lesion growth, and coalescence. These processes are supposed to capture the essence of TB dynamics in lungs, but of course they are not the only ones [6]. In fact, no model could be complete [37]. This also allows the use of a fewer number of parameters, when compared with other systems biology approaches to the same problem [1], and thus provides more robustness to the fitting. The principal novelty of this model is the implementation of the bubble model in an explicit space like the bronchial tree in order to simulate the endogenous reinfection processes. Nevertheless, it still has a few limitations that should be mentioned, the most important being the following: - Exogenous reinfection is not yet considered in this model. Its incorporation may change the outcome when simulating an ATB infection, as it acts as a new mechanism to generate new infection focuses. Nevertheless, the experimental data used in this study were obtained under conditions that prevented exogenous reinfection. Therefore, the inclusion of this mechanism should be supported by experimental designs that allow it. - The bronchial tree model is absolutely deterministic, for now. In the future we expect to add some random noise in this algorithm in order to obtain different bronchial trees from a single pulmonary surface. This will be useful for analysing the role of specific bronchial tree properties in TB evolution as well as to account for heterogeneity sources. - Infection spreading parameters are uniform in each virtual lung. Nevertheless, breathing amplitude is not constant, but varies from lower and middle lobes (wider amplitude) to upper ones (lower amplitude). Breathing amplitude is probably related with lesion spreading; then, higher values of βwould better fit the local behaviour of less dispersion in the upper lobe, while lower values of βwould be appropriate for describing the spreading in the lower and PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 20 / 25 middle lobes. This should be taken into account by creating a βvariable profile inside the bronchial tree. With our model, we have seen that a bigger dispersion reduces the probability of causing big lesions, which is true locally. At the same time, greater dispersion may increase the probability of generating new infection foci in parts of the lung with a smaller breathing amplitude, and therefore with a higher probability of evolving towards bigger granulomas. These dynamics will be carefully explored in the future. Finally, we must mention the four principal assumptions that could be refined and even refuted in the future: a. Our model follows dynamic hypothesis assumptions [38], but there are more hypotheses that can be submitted to check feasibility such as [39,40]. In addition, other important processes such as the role of oxygenation [41,42] could be incorporated. b. The model assumes that the larger the mother, the more likely it is to generate new lesions. This is based on the assumption that bigger lesions would have more foamy macrophages that are more likely to be drained and, therefore, to be able to cause a reinfection [43]. It is also supported by experimental observations in macaques of a higher tendency of bigger lesions to disseminate [44]. c. The reinfection model also considers that the older a lesion is, the less likely an infected macrophage is to escapes from it. This is in accordance with experimental observations [22]. In the future, a term that depends on the initial infection time in addition to the age of the lesion could be taken into consideration because, as can be seen in [45], initial lesions are bigger than their daughters in ATB. d. An infection is considered to correspond to an active disease if a lesion grows larger than 1 cm. ATB definition is more complicated than simply having a lesion bigger than 1 cm [35]. Then, other indicators should be explored to define an ATB. In fact, this threshold is not a standard accepted boundary. In supplementary material, section 5 (Fig E in S1 File), an exploration of different values of this threshold is shown. The strength of this model is that the same tendencies that were observed in Fig 10 can be observed in the different exploration threshold simulations. In general, model falsifiability could be successfully carried out with data from CT (or equivalent) of TB dynamics in big mammals, at different time points. The fact that only a final photo of the system is available is clearly a drawback for the testing of the processes involved. Indeed, tracking the growth of individual lesions at several timepoints should be enough for testing the generalized logistic model. Barcoding techniques have also shown their appropriateness to distinguish contained from disseminated lesions in macaques [44], and thus could provide a way to determine which the initial granulomas are and a pattern of dissemination. Bacteria barcoding, together with a time tracking for the 3D characterization of location and size of lesions by means of CT, can provide key information about the range and relative importance of dissemination, as well as the possible geometrical constraints. Nevertheless, one of the drawbacks of the required experimental tests that should be always kept in mind is that the n is usually small, while the intraspecific diversity is high. This approach has consisted of the testing of a single model, which seems to go against the strong inference in mathematical modelling [46]. Nevertheless, we have focused on exploring and exploiting all the possibilities given by this model and the available experimental data. In addition to the above-mentioned search for new experimental measurements that can refute the stated hypotheses, future work should include the testing of alternative models whose rejection would provide more clues on the natural history of tuberculosis. PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 21 / 25 Conclusions We built a computational model which includes the virtual lung and the updated Bubble model. This model successfully fits CT experimental observations of latent tuberculosis in minipigs. The model incorporates the basis of the dynamic hypothesis [7,37] to reproduce lesion propagation through the bronchial tree on a latent tuberculosis infection. The agreement between experimental and numerical results reinforces the feasibility of the dynamic hypothesis, i.e., it is able to explain the experimental results observed in latently infected minipigs. In particular, we have observed the importance of the bronchial tree in the endogenous reinfection (β>0). Most important parameters of the model could be related with the corresponding biophysical processes. Therefore, the model is consistent with the data. Parameter βcan be related with the breathing amplitude as a factor determining how far new lesions can appear; r max may be related with the effect of inflammatory response of the host, as it is the main cause of lesion growth; and ρdetermines the probability of triggering the endogenous reinfection process. It has been seen that ν, the growth velocity of the lesions, can play a similar role in triggering endogenous reinfection process, while slowing the growth velocity of lesions may be a mechanism for detaining endogenous reinfection. According to the in silico experiments carried out with the model, an active TB in an immunocompetent host may be caused by high inflammatory response or by moderate inflammatory response combined with small breathing amplitude. The former has already been observed experimentally in several studies [15,21,22]. The latter may be a clue for understanding the usual presence of active TB in the upper lobe, as suggested by Cardona & Prats [9], since it is in the lung zone that the breathing amplitude is smaller. In fact, further refinements and updates of the model should include inter-lobular differences so that this last possibility can be carefully explored. Supporting information S1 Dataset. Experimental dataset of CT measurements. Sheet “Lungs” contains, for each minipig, the pulmonary geometrical measurements and the global distribution of lesions in each lung. Sheet “Lesions” contains the measured characteristics of each lesion. (XLSX) S1 File. Supplementary material. This file includes the following supplementary sections: (1) Sensitivity analysis for transition and active set; (2) Extended sensitivity analysis; (3) Diameter-length relation in bronchial tree; (4) Details of the simplified protocol for model parametrization sensitivity; and (5) Thresholds exploration. (PDF) Acknowledgments We sincerely thank Marina Vegue ´and Carles Barril for the building of the preliminary human bronchial trees, which was the basis for the current algorithm. We thank Bernat Puig for his help on building the computational pulmonary surface. Author Contributions Conceptualization: Jordi Bechini, Montserrat Tenesa, Ricardo Pe ´rez, Cristina Vilaplana, Joaquim Valls, Daniel Lo ´pez, Pere-Joan Cardona, Clara Prats. PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 22 / 25 Data curation: Martı ´Català, Jordi Bechini, Montserrat Tenesa, Ricardo Pe ´rez, Mariano Moya, Pere-Joan Cardona. Formal analysis: Martı ´Català, Sergio Alonso, Daniel Lo ´pez, Clara Prats. Funding acquisition: Pere-Joan Cardona, Clara Prats. Investigation: Martı ´Català, Jordi Bechini, Montserrat Tenesa, Ricardo Pe ´rez, Mariano Moya, Cristina Vilaplana, Pere-Joan Cardona, Clara Prats. Methodology: Martı ´Català, Jordi Bechini, Montserrat Tenesa, Ricardo Pe ´rez, Cristina Vilaplana, Daniel Lo ´pez, Pere-Joan Cardona, Clara Prats. Resources: Cristina Vilaplana, Pere-Joan Cardona, Clara Prats. Software: Martı ´Català, Joaquim Valls, Daniel Lo ´pez, Pere-Joan Cardona, Clara Prats. Supervision: Cristina Vilaplana, Pere-Joan Cardona, Clara Prats. Validation: Martı ´Català, Pere-Joan Cardona, Clara Prats. Visualization: Martı ´Català, Joaquim Valls, Sergio Alonso, Clara Prats. Writing – original draft: Martı ´Català, Daniel Lo ´pez, Pere-Joan Cardona, Clara Prats. Writing – review & editing: Martı ´Català, Jordi Bechini, Montserrat Tenesa, Ricardo Pe ´rez, Cristina Vilaplana, Joaquim Valls, Sergio Alonso, Daniel Lo ´pez, Pere-Joan Cardona, Clara Prats. References 1. 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 2. World Health Organization. Global Tuberculosis report. Technical report. 2018. Available from: https:// www.who.int/tb/publications/global_report/en/ 3. Behr MA, Edelstein PH, Ramakrishnan L. Revisiting the timetable of tuberculosis. BMJ. 2018 Aug 23; 362:k2738. https://doi.org/10.1136/bmj.k2738 PMID: 30139910 4. Bermudez LE, Danelishvili L, Early J. Mycobacteria and macrophage apoptosis: complex struggle for survival. Microbe. 2006; 1(8):372–375. https://doi.org/10.1128/microbe.1.372.1 5. Mahamed D, Boulle M, Ganga Y, Mc Arthur C, Skroch S, Oom L, et al. Intracellular growth of Mycobacterium tuberculosis after macrophage cell death leads to serial killing of host cells. Elife. 2017; 6: e22028. https://doi.org/10.7554/eLife.22028 PMID: 28130921 6. Cardona PJ. Pathogenesis of tuberculosis and other mycobacteriosis. Enfermedades Infecciosas y Microbiologı ´a Clı ´nica. 2017; 36(1):38–46. https://doi.org/10.1016/j.eimc.2017.10.015 PMID: 29198784 7. Cardona PJ. Revisiting the natural history of tuberculosis. The inclusion of constant reinfection, host tolerance, and damage-response frameworks leads to a better understanding of latent infection and its evolution towards active disease. Arch Immunol Ther Exp (Warsz). 2010; 58(1): 7–14. https://doi.org/ 10.1007/s00005-009-0062-5 PMID: 20049645 8. Verrall A, Netea M, Alisjahbana B, Hill P, van Crevel R. Early clearance of Mycobacterium tuberculosis: a new frontier in prevention. Immunology. 2014 Apr; 141(4):506–13. https://doi.org/10.1111/imm.12223 PMID: 24754048 9. Cardona PJ, Prats C. The small breathing amplitude at the upper lobes favors the attraction of polymorphonuclear neutrophils to Mycobacterium tuberculosis lesions and helps to understand the evolution toward active disease in an individual-based model. Front Microbiol. 2016 Mar 29; 7:354. https://doi. org/10.3389/fmicb.2016.00354 PMID: 27065951 10. 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 PMID: 26986986 11. Ibargu¨en-Mondrago ´n E, Esteva L, Burbano-Rosero EM. Mathematical model for the growth of Mycobacterium tuberculosis in the granuloma. Math Biosci Eng. 2018 Apr 1; 15(2):407–428. https://doi.org/ 10.3934/mbe.2018018 PMID: 29161842 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 23 / 25 12. Segovia-Juarez JL, Ganguli S, Kirschner D. Identifying control mechanisms of granuloma formation during M.tuberculosis infection using an agent-based model. J Theor Biol. 2004 Dec 7; 231(3):357–76. https://doi.org/10.1016/j.jtbi.2004.06.031 PMID: 15501468 13. Marino S, Kirschner D. A multi-compartment hybrid computational model predicts key roles for dendritic cells in tuberculosis infection. Computation (Basel). 2016; 4(4). pii: 39. https://doi.org/10.3390/ computation4040039 PMID: 28989808 14. Prats C, Vilaplana C, Valls J, Marzo E, Cardona PJ, Lo ´pez D. Local inflammation, dissemination and coalescence of lesions are key for the progression toward active tuberculosis: the Bubble model. Front Microbiol. 2016 Feb 2; 7:33. https://doi.org/10.3389/fmicb.2016.00033 PMID: 26870005 15. 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 (Edinb). 2014 Jan; 94(1):55–64. https://doi. org/10.1016/j.tube.2013.09.004 PMID: 24291066 16. Richards F. A flexible growth function for empirical use. Journal of experimental Botany. 1959; 10 (2):290–301. 17. Peake JL, Pinkerton KE. Gross and subgross anatomy of lungs, pleura, connective tissue septa, distal airways, and structural units. In: Comparative biology of the normal lung 2015 Jan 1 (pp. 21–31). Academic Press. 18. Cardona PJ. The key role of exudative lesions and their encapsulation: lessons learned from the pathology of human pulmonary tuberculosis. Front Microbiol. 2015 Jun 16; 6:612. https://doi.org/10.3389/ fmicb.2015.00612 PMID: 26136741 19. North RJ, Jung YJ. Immunity to tuberculosis. Annu. Rev. Immunol. 2004 Apr 23; 22:599–623. https:// doi.org/10.1146/annurev.immunol.22.012703.104635 PMID: 15032590 20. Bailey M, Christoforidou Z, Lewis MC. The evolutionary basis for differences between the immune systems of man, mouse, pig and ruminants. Veterinary immunology and immunopathology. 2013 Mar 15; 152(1–2):13–9. https://doi.org/10.1016/j.vetimm.2012.09.022 PMID: 23078904 21. Ramos L, Obregon-Henao A, Henao-Tamayo M, Bowen R, Lunney JK, Gonzalez-Juarrero M. The minipig as an animal model to study Mycobacterium tuberculosis infection and natural transmission. Tuberculosis (Edinb). 2017 Sep; 106:91–98. https://doi.org/10.1016/j.tube.2017.07.003 PMID: 28802411 22. Gil O, Diaz I, Vilaplana C, Tapia G, Diaz J, Fort M, et al. Granuloma encapsulation is a key factor for containing tuberculosis infection in minipigs. PLoS One. 2010 Apr 6; 5(4):e10030. https://doi.org/10. 1371/journal.pone.0010030 PMID: 20386605 23. Bechini J. Estudio de la tuberculosis pulmonar mediante Tomografı ´a Computarizada Multidetector en un modelo experimental de minipig. PhD thesis. Universitat Autònoma de Barcelona. 2016. Available from: https://ddd.uab.cat/pub/tesis/2017/hdl_10803_400765/jbb1de1.pdf 24. CatalàM. Modelling and simulation of tuberculosis lesions dynamics in a minipig bronchial tree. Bachelor thesis. Universitat Politècnica de Catalunya. 2015. Available from: https://upcommons.upc.edu/ handle/2117/77133 25. Vegue ´M. Model tridimensional de l’arbre bronquial humàper a l’estudi de la dissemeniacio ´de Mycobacterium tuberculosis. Master’s thesis. Universitat de Barcelona. 2012. 26. Weibel E. Morphometry of the human lung. New York (USA). Academic Press; 1963. 27. Singh VK, Thrall KD, Hauer-Jensen M. Minipigs as models in drug discovery. Expert Opin Drug Discov. 2016 Dec; 11(12):1131–1134. https://doi.org/10.1080/17460441.2016.1223039 PMID: 27546211 28. Murray CD. The physiological principle of minimum work: I. The vascular system and the cost of blood volume. Proc Natl Acad Sci USA. 1926 Mar; 12(3):207–14. https://doi.org/10.1073/pnas.12.3.207 PMID: 16576980 29. Florens M, Sapoval B, Filoche M. An anatomical and functional model of the human tracheobronchial tree. J Appl Physiol (1985). 2011 Mar; 110(3):756–63. https://doi.org/10.1152/japplphysiol.00984.2010 PMID: 21183626 30. Rozanek M, Roubik K. Mathematical model of the respiratory system—comparison of the total lung impedance in the adult and neonatal lung. International Journal of Biomedical Sciences 2007; 249–252. 31. Murray CD. The physiological principle of minimum work applied to the angle of branching arteries. J Gen Physiol. 1926 Jul 20; 9(6): 835–841. https://doi.org/10.1085/jgp.9.6.835 PMID: 19872299 32. Kitaoka H, Takaki R, Suki B. A three-dimensional model of the human airway tree. J Appl Physiol (1985). 1999 Dec; 87(6):2207–17. 33. Ginovart M, Prats C, Portell X. Microbial individual-based models and sensitivity analyses: local and global methods. In: International Conference on Predictive Modeling in Foods (Dublin), 2010; pp. 313– 316. Available from: https://upcommons.upc.edu/handle/2117/13546 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 24 / 25 34. Marino S, Hogue IB, Ray CJ, Kirschner DE. A methodology for performing global uncertainty and sensitivity analysis in systems biology. Journal of theoretical biology. 2008 Sep 7; 254(1):178–96. https://doi. org/10.1016/j.jtbi.2008.04.011 PMID: 18572196 35. Andreu J, Ca ´ceres J, Pallisa E, Martinez-Rodriguez M. Radiological manifestations of pulmonary tuberculosis. Eur J Radiol, 51 (2004), pp. 139–149 https://doi.org/10.1016/j.ejrad.2004.03.009 PMID: 15246519 36. McCullagh P, Nelder JA. Generalized linear models, 2nd edn. ( Chapman and Hall: London). Standard book on generalized linear models. 1989. 37. Oreskes N, Shrader-Frechette K, Belitz K. Verification, validation, and confirmation of numerical models in the earth sciences. Science. 1994 Feb 4; 263(5147):641–6. https://doi.org/10.1126/science.263. 5147.641 PMID: 17747657 38. Cardona PJ. A dynamic reinfection hypothesis of latent tuberculosis infection. Infection. 2009 Apr 1; 37 (2):80. https://doi.org/10.1007/s15010-008-8087-y PMID: 19308318 39. Gong C, Linderman JJ, Kirschner D. A population model capturing dynamics of tuberculosis granulomas predicts host infection outcomes. Mathematical biosciences and engineering: MBE. 2015 Jun; 12 (3):625. https://doi.org/10.3934/mbe.2015.12.625 PMID: 25811559 40. Cohen SB, Gern BH, Delahaye JL, Adams KN, Plumlee CR, Winkler JK, et al. Alveolar macrophages provide an early Mycobacterium tuberculosis niche and initiate dissemination. Cell host & microbe. 2018 Sep 12; 24(3):439–46. https://doi.org/10.1016/j.chom.2018.08.001 PMID: 30146391 41. Sershen CL, Plimpton SJ, May EE. Oxygen modulates the effectiveness of granuloma mediated host response to Mycobacterium tuberculosis: a multiscale computational biology approach. Frontiers in cellular and infection microbiology. 2016 Feb 15; 6:6. https://doi.org/10.3389/fcimb.2016.00006 PMID: 26913242 42. Bowness R, Chaplain MA, Powathil GG, Gillespie SH. Modelling the effects of bacterial cell state and spatial location on tuberculosis treatment: Insights from a hybrid multiscale cellular automaton model. Journal of theoretical biology. 2018 Jun 7; 446:87–100. https://doi.org/10.1016/j.jtbi.2018.03.006 PMID: 29524441 43. Russell DG, Cardona PJ, Kim MJ, Allain S, Altare F. Foamy macrophages and the progression of the human tuberculosis granuloma. Nature immunology. 2009 Sep; 10(9):943. https://doi.org/10.1038/ni. 1781 PMID: 19692995 44. Martin CJ, Cadena AM, Leung VW, Lin PL, Maiello P, Hicks N, et al. Digitally barcoding Mycobacterium tuberculosis reveals in vivo infection dynamics in the macaque model of tuberculosis. MBio. 2017 Jul 5; 8(3):e00312–17. https://doi.org/10.1128/mBio.00312-17 PMID: 28487426 45. Sharpe S, White A, Gleeson F, McIntyre A, Smyth D, Clark S, et al. Ultra low dose aerosol challenge with Mycobacterium tuberculosis leads to divergent outcomes in rhesus and cynomolgus macaques. Tuberculosis. 2016 Jan 1; 96:1–2. https://doi.org/10.1016/j.tube.2015.10.004 PMID: 26786648 46. Ganusov VV. Strong Inference in Mathematical Modeling: A Method for Robust Science in the TwentyFirst Century. Front. Microbiol. 2016; 7:1131. https://doi.org/10.3389/fmicb.2016.01131 PMID: 27499750 PLOS COMPUTATIONAL BIOLOGY Dynamics of TB lesions in a virtual lung PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007772 May 20, 2020 25 / 25