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

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 agent-based 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.


Sensitivity analysis for transition and active set
The sensitivity analysis is performed using as initial infection a single lesion on the mass centre of the experimental lesions observed experimentally in each minipig.
The sensitivity analysis is performed for 5 variables: , , , , and . A set of 11 simulations with different parameter values is designed: one simulation with all the parameters at their default values, five simulations increasing one-at-a-time parameter by 10%, and five simulations decreasing one-at-a-time parameter by 10%. We compute an ANOVA test analysis for 4 outcomes by comparing the results obtained with the original set of parameters and those obtained with each new parameter combination [1].
In the main article, sensitivity analysis is performed for latent set of parameters: ρ = 0.12 day -1 , β = 0.08 mm -1 , and rmax = 0.68 mm (Table 4). In Table A, the results for a transition set of parameters are shown: ρ = 0.045 day -1 , β = 0.10 mm -1 , and rmax = 6.5 mm. In Table B, the results for an active set of parameters are shown: ρ = 0.045 day -1 , β = 0.15 mm -1 , and rmax = 10 mm. These tables show that sensitivity analysis depends on the point chosen, but there are some relations that can be observed for the three sets used. There is a dependence between the number of lesions and simulation time (Tmax). There is also a dependence between mean diameter and β, rmax, ν, and Tmax; dispersion and β show a dependence, as well. Finally, a dependence between the number of coalescences and ρ and ν is found, but no relation with Tmax.
Other relations are significant in just one or two of the points where the analysis is carried on. For example, in transition and active sets a relation between the number of lesions and β and rmax is found, which is not observed in the latent set.
A point that may be noted is that sensitivity for ρ and ν parameters is almost the same in the three analysed points. Thus, in this model ρ and ν have a very similar role. 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.

Extended sensitivity analysis
An uncertainty and sensitivity analysis is performed for our model as described in [2] using a sampling-based method. Thus, 1000 different parameter sets are used to explore the space using a Latin Hypercube Sampling (LHS) technique.
Parameters are explored in the ranges: Where TH is the threshold that determines if an infection is considered an Active Disease indicator increases with rmax, ρ, ν, Tmax, and β, while it decreases with TH. In fact, as can be seen in Fig A, the only significant relation of TH is with disease indicator outcome. This is a coherent result because TH is the threshold that determines whether an infection is considered ATB or not.
Maximum diameter relations are the same as disease indicator relations. This is a coherent result because the disease indicator is directly related with the maximum diameter of the observable lesions.
The parameter that most affects the dispersion is β. Larger values of β entail a decrease in dispersion, as can be seen in Fig 4A of  Coalescence process is affected by all parameters except TH, as mentioned before.
In fact, it is a crucial process that is essential in explaining most of the relations seen in this sensitivity analysis. The two parameters that most affect the coalescence number are ν and ρ due to their relation with the number of lesions.

Diameter-length relation in bronchial tree
Although length and diameter of branches in a bronchial tree have a complex relation, a common approximation used in computational bronchial trees is that the length of a certain branch is 3 times its diameter [3,4].
We used Rozanek and Roubik's [5] data from an adult human bronchial tree to determine the linear relation between diameter and length of the different branches. Using linear regression as shown in Fig B, length and diameter relation can be approximated by: With a goodness of fit: This reinforces the validity of the approximation used in our model.

Details of the simplified protocol for model parametrization sensitivity
We aimed to design a protocol for parameterization in order to find the set of parameters (rmax, ρ and β) that best fit experimental data. We started our simulations with a single lesion at the mass centre of the experimentally observed distribution. To fit the three remaining parameters, we built three error functions to evaluate the agreement of three outcomes of the computer simulations with experimental data: NLE, DE, and SE, as explained in the main article.
In order to determine how each input parameter affects the defined errors, we ran These cases correspond to the green cells in Table 3

Threshold exploration
ATB is determined by the biggest lesion. As seen in [6], lesions bigger than 1 cm can be observed in a thorax radiography. Then, threshold radius (TH) is fixed at 5 mm. In Fig 8 of