Application of a Stochastic Modeling to Assess the Evolution of Tuberculous and Non-Tuberculous Mycobacterial Infection in Patients Treated with Tumor Necrosis Factor Inhibitors

In this manuscript we apply stochastic modeling to investigate the risk of reactivation of latent mycobacterial infections in patients undergoing treatment with tumor necrosis factor inhibitors. First, we review the perspective proposed by one of the authors in a previous work and which consists in predicting the occurrence of reactivation of latent tuberculosis infection or newly acquired tuberculosis during treatment; this is based on variational procedures on a simple set of parameters (e.g. rate of reactivation of a latent infection). Then, we develop a full analytical study of this approach through a Markov chain analysis and we find an exact solution for the temporal evolution of the number of cases of tuberculosis infection (re)activation. The analytical solution is compared with Monte Carlo simulations and with experimental data, showing overall excellent agreement. The generality of this theoretical framework allows to investigate also the case of non-tuberculous mycobacteria infections; in particular, we show that reactivation in that context plays a minor role. This may suggest that, while the screening for tuberculous is necessary prior to initiating biologics, when considering non-tuberculous mycobacteria only a watchful monitoring during the treatment is recommended. The framework outlined in this paper is quite general and could be extremely promising in further researches on drug-related adverse events.


Introduction
Over the last decades the improved understanding of the pathogenesis of chronic inflammatory diseases, together with a major advance in biotechnology, have accelerated the development of biological therapies, designed to neutralize specific targets that mediate and sustain the clinical manifestations of diseases. These compounds, mainly monoclonal antibodies (mAb) and fusion proteins, introduced a breakthrough in the management of different conditions including inflammatory rheumatologic disorders [1]. In this context, the first setting of application of the biological agents was rheumatoid arthritis (RA), a chronic autoimmune disease affecting approximately 1% of the adult population [2]. If the disease is not treated adequately, progressive deformity can lead to loss of quality of life and reduce average life expectancy by about a decade [2]. Studies on the pathogenic mechanisms of RA have revealed that tumor necrosis factor (TNF) is a cytokine playing a critical role in the inflammatory cascade that results in the irreversible joint damage typical of the disease [3]. Following these discoveries, a series of clinical trials in patients with RA showed the therapeutic benefit of TNF blockade [4]. As a consequence, five biological agents engineered to block TNF actions are currently available: infliximab, adalimumab, golimumab, certolizumab pegol (all of them mAb), and etanercept (a receptor fusion protein) [5]. While being highly effective, TNF blockers have raised concerns about the potential for an increased susceptibility to infections, in particular the reactivation of latent tuberculosis (TB) infection [6][7][8][9][10]. Mycobacterium tuberculosis, the cause of human TB, can result in a metastable clinical latency lasting for decades. Much has been speculated about the structure of granuloma which should contain Mycobacteria, since murine models indicated that TNF was necessary for both formation and maintenance of granulomas [11]. However, subsequent studies on zebrafish model [12], monkeys [13], and humans [14][15][16][17] challenged these data, demonstrating that the crucial role of TNF in the granuloma was indeed macrophage activation and stimulation of chemokine production. The reactivation of latent TB infection has been associated with all TNF inhibitors, hence pre-initiation screening procedures have been recommended, which have successfully reduced the number of reported cases [18], although current screening tools lack sensitivity and specificity [19,20].
TNF blockers seem to increase also the risk of other granulomatous diseases, but little is known about the emergence of illness due to non-tuberculous Mycobacteria (NTM). These are a huge ensemble of pathogens (e.g. M. avium, M. abscessus, and so on [21]) and up to date, approximately 50 different mycobacterial species are considered to be etiological agents of human diseases and this number seems still growing [21]. Most cases these days occur in hosts with relatively intact immune responses. However, RA and other chronic diseases with pulmonary manifestations can predispose a person to NTM pulmonary disease [22] expressing as a possible serious complication, especially in immunosuppressed subjects. Thus, it is of utmost importance to study also the risk related to NTM, in the perspective to understand if a proper screening may be helpful in conferring a wider protection to the patients. This is particularly true, in that the TNF blockers appear to predispose both to disseminated and localized disease [21,23], but also because these infections are increasing in prevalence, especially among women, which are more frequently affected by RA than men. In the present state of the art, the scenario for NTM diseases, with respect to TNF-blocking drugs, seems different from the TB counterpart: in particular, through extensive experimental screening, both Wallis and coworkers [10], and Winthrop and coworkers [23,24] evidenced that latency and reactivation do not seem to play a crucial role in this context, yet a clear-cut picture is still missing. Now, as far as TB is concerned, data collected through the Adverse Event Reporting System of the US Food and Drug Administration (FDA) in the time-window 1998{2002, related to the two test-case drugs with a different mechanism of action, i.e. infliximab and etanercept, highlight that TB infections involve 54 over 10 5 patients treated with infliximab and 28 over 10 5 patients receiving etanercept [10]. Therefore, the question is: As the latency in TB can last decades, are these infections (in patients under therapy) new ones or are they reactivation of previously encountered pathogens due to a suppressed immune system? This kind of question underlies the awareness of a real need and disposal for extensive pre-screening procedures. Unfortunately, the answer is by far not trivial as, for TB, there are no secure pathways to discriminate between a new infection or the raise of a previous one. Moreover, a clear methodology for finding latencies is still lacking. Furthermore, the rarity and different sizes of this infection in different countries (ranging from 5 over 10 5 in Sweden up to 140 over 10 5 in Romania [8]) implies that data analysis and its subsequent interpretation must be carefully performed.
As for NTM, still from FDA, through the post-marketing surveillance system (MedWatch) and through a further survey within the Emerging Infections Network of the Infectious Diseases Society of America (IDSA), Winthrop and coworkers reported a detailed study of possible correlations between the usage of TNF blockers and the emergence of NTM diseases: over a time-window of 8 years, they highlighted a higher prevalence of NTM diseases in patients treated with infliximab rather than etanercept [24].
In order to investigate possible correlations between the incidence of infections by such Mycobacteria and biological therapy, one could rely directly on the molecular details of TNF processing signal (which has been, at least partially, elucidated, see e.g. [25]), coupled to the underlying infliximab and etanercept mechanism of action, which could be achievable directly through molecular immunology approach. Beyond these ''standard'' strands, a completely different route can also be performed: Given the relative large amount of collected data, the problem can be considered from a purely inferential viewpoint, by-passing the underlying molecular immunology know-how (see also [26,27]). According to this perspective, in Ref. [9,10], an abstract (logical) environment for TB case has been defined, where patients can occupy one of the (following) five different states: (0) No infection, (1) New infection, (2) Latency, (3) Reactivation of a previous TB infection, (4) Post first TB encounter. Clearly, the patients starting the therapy (and hence belonging to the survey) can correspond to either state (0) or (2), because all the other states imply quantifiable sickness and the patient would then be treated for TB rather than RA. Then, at the end of the survey, a fraction of these patients will be in an illness state, i.e. either state (3) or (4). The transition rates between different states are assumed as freeparameters, whose values are estimated through numerical simulations: the best estimate is the one able to reproduce, with the smallest error, the experimental data. Remarkably, the probability of latent TB reactivation in patients treated with infliximab turned out to be an order of magnitude per unit of time higher than the same probability for patients trated with etanercept [9].
Here, we first formalize this approach in terms of Markov chains and we write the related Master equation in continuoustime limit, then we solve the model analytically and study its properties in full details. In this way we get the explicit expression for the number of patients c(t) exhibiting a TB (re)activation, as a function of time t. One step forward, we check the robustness of our results through extensive Monte Carlo simulations and over the clinical data of the TB scenario, finding overall excellent agreement among all our results (and previous literature). Moreover, we find that different magnitudes for the probability of reactivation correspond to qualitative different behaviors for c(t) (on the proper timescale), that is, the number of patients displaying active infection increases exponentially in time when using infliximab and linearly in time when using etanercept.
The analytical expression for the whole evolution of the system implies a great feasibility of the technique itself (e.g. we have the whole set of first integrals and a clear picture of all the hidden symmetries) and also allows to address, in complete generality, several instances. In particular, we can finally consider generic NTM infections, where, interestingly, the scenario appears quite different from the TB counterpart: clinical data suggest that c(t) (on the proper timescale) grows quadratically with time and this is recovered by our analytical picture only under the assumption of a negligible role played by latency reactivation. We check these findings also through extensive Monte Carlo runs, which are in full agreement too. Remarkably, this is very consistent with the present state of the art in the medical literature dealing with NTM.
As a final result, there are two types of conclusions which stem from our work: The former belongs to the world of modelers, while the latter to the world of clinicians.
From a mathematical perspective, the encouraging results of this approach may pave the way for the development of handily and fruitful instruments for physicians.
Much more carefully, in the clinician's counterpart, as this approach bypasses the whole underlying biological complexity, it may contribute to confirm, from a theoretical perspective, the current understanding of adverse events coupled to TNFinhibitors and the consequent real need for screening procedures before undergoing biological therapies.

Materials and Methods
In this section we formalize the scheme introduced in [9] and aimed to reproduce data of TB onset in patients treated with TNF inhibitors, with particular attention on infliximab (an anti-TNF mAb) and etanercept (a soluble TNF receptor). Seeking for clarity, in this section we mention only applications to the TB case, although, as we will see in the second part of the paper, the approach is rather robust and can be successfully applied to the NTM case, too.
The model, whose structure is depicted in Figure 1, consists in identifying a set of possible states for the patient subjected to biological treatments, and in fixing the likelihood for the patient to change his/her state within a proper unit time. Of course, on large samples, some patients may experience sudden incidents (e.g. death for other causes) or some others may assume both the drugs: the analysis has been previously purified from these cases [9,10].
The clinical states available to a test-patient are taken as follows (see Figure 1): 0 : Absence of infection; 1 : New infection (that after a time t can give rise either to active TB or latent infection); 2 : Latent infection; 3 : Reactivated TB after latency; 4 : Active TB (that progresses from new infection within a time t, without an intervening period of latency).
Moreover, each patient is assumed to change his/her state, following the corresponding transition probabilities, which constitute the model parameter set, and are meant over a proper unit time t. Using t to label the time, these probabilities are: L : Probability of having a latent infection at the beginning of the observation (t~0), while, obviously, (1{L) is the probability of not having any infection at that moment; N : Probability of TB infection during the observational time; P : Probability of a new TB infection to become active TB after a time t; as a consequence, (1{P) is the probability of this new infection to give rise to a latent infection; R : Probability of reactivation of a latent TB infection.
We stress that such probabilistic framework is based on purely clinical variables.
On the experimental side, the available data consist in a collection of times (one for each patient) corresponding to the onset of TB (in its active phase, namely a detectable scenario), after the beginning of the treatment with TNF blockers. As a consequence, the only states which are possible to observe are the states 3 and 4. Unfortunately, as discussed in the introduction, these states (that account for ill patients) are not distinguishable one respect to the other by simply looking at the data (hence motivating both earlier studies [9,10] and our machinery), however, some progress can be made using stochastic extremization. The idea resembles the standard maximum likelihood and consists in finding the best values for free parameters such that the theoretical curves collapse over the experimental data [28].
As already outlined, following this procedure, the main result in [9] is that the principal difference between infliximab and etanercept treatments resides in different management of latent TB, i.e. on R. The former seems to enhance reactivation one order of magnitude more than the latter.

Markov Chains and Master Equations
The model described in the previous section can be translated into a set of differential equations coding for the temporal evolution of the probability of patient's states (which can be compared to the corresponding fractions over a sample of patients given the large collection of data).
Being the states discrete, this can be accomplished in complete generality using Markov chains, namely a (discrete-time) probabilistic framework where the probability of being in a given state at a given time t depends only on the probability distribution over all the states at the previous time step t{1, and on the transition rates linking these states.
It is instructive to consider the illustrative Markov chain with only three states (A, B and C), non-null transition rates w A?B and w B?C and time step Dt, shown in Figure 2.
Note that, in the model, the probabilities of going from A to B and from B to C exist but not the opposite (w B?A~0 , w C?B~0 ) hence, if the initial state is all concentrated in C, there will be no evolution, while if the starting point is spread among A and B, after enough time, the probability distribution will be peaked on C only (but in its finite temporal evolution resides our interest). Now, the probability of remaining-at/moving-into the state B in the time interval Dt is given by the probability of already being in B (hence p B (t)) plus the probability of arriving in B from A times the probability of being in A at the previous step (hence w A?B p A (t)) minus the probability of leaving B to C times the probability of being in B at the previous step (hence w B?C p B (t)); this concept can be written as follows: Now, the mathematics for continuous variable differential equations is much more handily and does not change significantly the perspective if the time step is small with respect to the global time window; indeed, both experimental data sets considered here (for TB cases [10] and for NTM cases [23]) fulfill this requirement and thus we are allowed to consider the time as a continuous variable. This can be achieved straightforwardly starting from the previous equation using a limit procedure: The evolution for the probability p B (t) is then ruled by the following differential equation, namely a ''Master equation'',  which acts as a continuous counterpart of the Markov chain in the discrete-time case: In general, for a system which can be in one of M arbitrary states, we need a M|M transition-rate matrix w (where w i?j is the rate for the transition from state i to state j) and the Master equation takes the form Note that we use the symbol : above the functions meaning their temporal derivative. Finally, we switch to a form where the explicit timescale t of the process appears directly in the equation, that is where W i?j *w i?j t stands for the probability of transition from state i to j along the time interval t.

Master Equations for the Model
Keeping in mind Figure 1, we can write down the system of differential equations describing the evolution of the five states earlier introduced as follows: with initial conditions p 0 (t~0)~1{L, p 2 (t~0)~L, p 1 (t~0)~p 3 (t~0)~p 4 (t~0)~0: The numbers indexing the probabilities mirror the enumeration of the previous section, that is, p 0 stands for the probability that a patient has never been affected by the infection, and so on. The parameter t represents the typical time for a patient experiencing a new infection to either develop the disease or to fall into a latent state and it should be chosen according to the natural time-scale of the process described. For instance, for the TB case, the data collected suggest that t is order of a few months [10,29], and we set t~1month for the sake of simplicity and in agreement with previous works [9,10]. Note that, as patients affected by active TB do not start RA therapy, we set p 1 (t~0)~p 3 (t~0)~p 4 (t~0)~0. Furthermore, the parameter L tunes the initial amount of latent-TB patients with respect to free-TB patients, such that for L~0 all patients are healthy, while for L~1 all patients display a latent TB infection; as we have no ways to discriminate between healthy and latentinfected patients, L is taken as a free parameter which can be estimated a posteriori comparing the solution of (5) with available data.
The solution of the system (5) can be easily obtained using first order ordinary differential equations theory and reads off as Of course, since the total amount of patients is conserved, C 0~p0 zp 1 zp 2 zp 3 zp 4 is an integral of motion, that is Beyond C 0 , the system ( ) admits another integral of motion C 1 , namely This means that the rate of growth for patients in the latency branch (i.e. in states 2,3) equals the rate of growth for the rest of infected patients (i.e. in state 4) weighted by a factor P {1 {1, so that the smaller P and the larger the difference between the related rates. The knowledge of integrals of motion can be very useful as they allow to obtain information in a very simple way; for instance, should P drop, then p 4 (t) would also decrease (or, analogously, p 2 (t)zp 3 (t) would increase) in order to maintain C 1 constant. Given C 0 and C 1 , other integrals of motion, which are combination of C 0 and C 1 , can be trivially built. For example, We underline that this kind of investigation can be accomplished only through an analytical study of the system. As discussed above, the fraction of active TB cases is given by the sum of the fraction of cases of direct TB after infection and of the fraction of cases with reactivated TB; namely, calling c(t) the total fraction of cases, we have: In the above quantity, the time dependence appears only through three different exponential decay terms (e {Nt=t , e {Rt=t , e {t=t ), which vanish at infinite time, so that the solution becomes a constant term equal to 1, meaning that, if we wait for a sufficient long (possibly infinite) time, all patients become sick (although, obviously, they can possibly die earlier due to reasons not related to RA/TB). In order to deepen the temporal evolution of these probabilities at relatively short times, it is useful to use a little bit of algebraic manipulation to distinguish constant terms from decaying terms, in such a way that we get where the three constants k 1 , k R and k N are related to the physiological parameters by Of course, from c(t) one can derive the effective number of cases multiplying c(t) by the overall number of treated patients. Before turning attention to the fitting procedure, we stress that the analytical solution in Eq. (7) was successfully checked through numerical methods, i.e. fourth-order Runge-Kutta algorithm and Monte Carlo simulations (notice that here, with Monte Carlo simulation we mean a simulation in which a set of virtual patients evolves in time following the Markov chain of Figure 1, giving a sample of the evolution of the fraction of cases during time. In our simulations we set 10 6 virtual patients).

The TB-infection Case
Having obtained the complete solution of the model and exploiting the available information on parameters, we now look for proper approximations able to highlight the effective behavior of c(t) in cases of practical interest, starting with the TB scenario.
In particular, Eq. (12) can be reduced to a simpler form if we note that the probability N of TB infection is much smaller than all the other parameters, in agreement with studies on TB and with results found in [9] and, a posteriori, in the current paper (see Table 1).
Hence, as a first approximation step, we assume N%1 and N=R%1 such that we can expand the solution, at the first order in N and N=R, as follows: notice that here and in the following we use the ''big-O'' Landau notation to characterize the growth rate of functions; more precisely, being f (x) and g(x) two arbitrary functions, we say f (x)~O(g(x)) as x?0 if there exists a positive real number M such that Df (x)DƒMDg(x)D.
By plugging the previous expressions into Eq. 12, with some algebra and retaining only up-to-linear terms in N or N=R, we get Let us now move further and focus on the exponential terms. First, we notice that 1wRwN and, consequently, we can neglect the term e {t=t , as it decays much faster that both e {Rt=t and e {Nt=t . Moreover, since the time range considered is &30 months and N is expected to be %t=t&10 {2 , we can expand e {Nt=t as e {Nt=t &1{Nt=t, and considering only the leading dependence on t, we get As for e {Rt=t , a similar approximation (e {Rt=t &1{Rt=t) can be adopted as long as R v * 10 {2 so to obtain the following linear approximation.
Notice that a smaller (larger) estimate for t would simply require a stricter (softer) condition on N and on R for the related linear expansions to hold (on the same time range); the model would not be affected and the parameters coupled with time, i.e. N,R, would be accordingly rescaled. As shown in Figure 3, the approximation (16) is rather good only for etanercept-treated patients, for which the best fit yields R~2:24 : 10 {2 . On the other hand, if we consider infliximab-treated patients, the approximation (16) does not fit data, while using (15) we get a good overlap with data and the best fit yields R~2:12 : 10 {1 , confirming that now Rt=t is no longer small over the time window. All best-fit coefficients are reported in Table 1, and the related errors are shown in Figure 4.
We can estimate how sensitive c(t) is with respect to the system parameters by deriving its analytic expression (see Eq. 12) with respect to N,P,R,L, respectively; in this way we get that, in the regime N%(1,P,R,L), the most relevant parameter affecting the behavior of c(t) is R. Another argument in favor of this claim is that, in the zero approximation of the solution (i.e. neglecting even terms O(N)), P does not appear at all.
In order to get further insight on the effect of infliximab and of etanercept on TB incidence, in Fig. 5 we plotted the model predictions for the percentage of patients having TB because of new infections (p 4 (t)) or reactivations (p 3 (t)).
To summarize, our results confirm that, in the present context, the most important difference between therapies based on infliximab or etanercept is that the former enhances TB reactivation more than the latter, in fact, we found R INF *10R ETA , in agreement with [9]. Such a discrepancy implies even a qualitatively different behavior of c(t) over the time-window considered: the number of infliximab-treated patients experiencing a TB infections grows exponentially in time, while for etanercept-treated patient the growth is linear.

The NTM-infection Case
The analytical expression for c(t), (see Eqs. 7 and 11) holds for a general environment schematizable as in Fig. 1; in the last paragraph we used the details of TB infections to implement convenient approximations. We now turn to NTM-infections in patients treated with infliximab and etanercept and look for proper approximations able to highlight the characteristic features of such case. The experimental data we refer to are those reported in [23] and collected over the period 1999{2006. Overall, there were 239 reports of NTM infections in patients who were receiving anti-TNF therapy. Most reports were for patients receiving infliximab (75%) or etanercept (17%) and here we shall focus just on these drugs.
Some remarks are in order here to merge mathematically the TB and NTM scenarios.
First, it is important to notice that patients affected by NTM lung infections typically suffer through long periods of illness before a clinical diagnosis is made. To comply with the actual state of the art on the involved time-scales we follow M. Iseman and T. Marras [29] that we quote: ''Preliminary prevalence estimates have been made, assuming that the disease duration for TB is 8 months and for pulmonary NTM is within the range of 4 to 10 years'' (48 to 120 months), hence for NTM the timescale is at least one order of magnitude larger than in the TB counterpart.
Another important point is that available data on NTM-diseases related to drug therapies lack the size of the survey, namely the number of patients participating to the screening, consequently, we can quantify the (cumulative) amount of sick patients, but we do not know their percentage. This is not a serious deficiency since, while we do not have access to the very reactivation probabilities R INF or R ETA , their ratio (which cancels out both time-scales and survey-size) still retains a quantitative information content.
As a result of the first remark, the (average) exit time from state 1, i.e. t, is comparable to the experimental time window, and one can expand the solution reported Eq. 7 at second order in t=t as We notice that, at this stage of expansion, p 4 (t) grows quadratically with time, while p 3 (t) presents two contributes growing, respectively, linearly and quadratically with time, the former being the leading one.
Such expressions must now be compared with experimental data, which display a purely quadratic growth in time, i.e. compatible with c(t)~at 2 (see Fig. 6). In order for the comparison to hold, one must therefore drop the linear coefficient, and this implies R*0 (or L*0 which is conceptually pretty similar). With this choice of parameters we get. and this form successfully fits experimental data, as shown in Figure 6; the best-fit coefficients are reported in Table 2, and the related errors are shown in Figure 7.
The functional form of p 4 (t) highlights why finding an estimate of the parameters is extremely difficult through blind numerical extremization, without any explicit analytical hint: The latencybranch turns out to be negligible in NTM-context and this results in a coupling of the parameters P and N as only a single probability streaming toward p 4 (t) survives. Therefore the two parameters alone are undetectable.
Using the coefficients of Table 2 in numerical simulations we obtain results in excellent agreement with real data and analytical ones. Notice that, impressively, the fit is very good although here the optimization relies only on one parameter.
We stress that the qualitative difference between the behaviors of c(t) in the case of TB and NTM infection (see Eqs. 15 and 19) mainly stems from a significant gap in the related disease duration (few months versus several years). In the latter case t is relatively large to allow the expansion of the exponential functions into a polynomial form, while for TB this can be accomplished only under the condition of small enough R. Such expansions do not significantly alter the predicted values of the fit parameters, as long as the underlying assumptions are consistent with facts.
To summarize, the theoretical framework developed evidences that the rate of reactivation R is vanishing: This issue makes the fraction p 3 of reactivation cases negligible, while the number of activation cases grows quadratically with time through the contribution of only p 4 , consistently with the actual understanding of NTM-pathology achieved through standard pathways. Moreover, differently from the TB case, in the NTM scenario no qualitative difference can be detected between infliximab and etanercept: The parabolic behavior for c(t) seems robust.

Discussion
In the last decades, several tools and concepts stemmed from the fields of stochastic mathematics and theoretical physics have been proposed to help the investigation of the immunological world, ranging from kinetic theories [30], to associative neural networks [31,32], to cellular automata and more [33,34].
Along the same line, in these notes we formalized and extended a stochastic approach to data analysis (originally introduced in [7-10]) for evidencing underlying correlations between adverse events and therapies based on immunosuppressants. In particular, the focus of our investigation concerns the risk of reactivation of latent mycobacterial infection in patients undergoing treatment with TNF-inhibitors.    We gave a clear and complete mathematical backbone to this approach, building it on explicit Markov processes, whose continuous-time limit yields the Master equation governing the evolution of the expected fraction of patients c(t) exhibiting an active infection. We also solved the Master equation in all details finding an analytical expression for c(t): Such mathematical developments make the original approach much more versatile and general: For instance, handling the complete (mathematical) solution allows to better account for reasonable approximations, tackling their control quantitatively (e.g. finding the proper timescales involved in the process or the integrals of motion constraining the evolution of the system). Furthermore, we can finally consider, within the same framework, different problems. In particular, we focused on TB and NTM infections emerging during anti-TNF therapies (infliximab and etanercept) according to data reported in [9,23].
In the former case, we recovered previous findings showing that the rate of reactivation R of TB from a latent state to an active state plays a crucial role: being R INF *10R ETA we get qualitatively different behaviors for c(t). More precisely, once fixed the observational time-window, for infliximab c(t) grows exponentially with time, while for etanercept it grows linearly with time. Hence, these results sustain the need, for patients candidate to TNF blockers, to perform an accurate TB screening at baseline, irrespective of the type of antiTNF. Indeed, screening may decrease the risk of TB reactivation in such patients, while it is less clear what should be done to prevent NTM disease occurrence or progression in patients taking biologic agents. Importantly, for this purpose, we found that the comparison with experimental data allows to infer that reactivation plays a very minor role for both the therapies and that c(t) grows quadratically with time.
We checked our results also against Monte Carlo simulation with excellent agreement.
Furthermore, our results are all consistent with recent experimental data and seem to indicate that TB and NTM infections are sustained by different pathogenetic mechanisms.
Non-tuberculous Mycobacteria are present in large numbers in the environment, including fresh water, aerosols, biofilms, and soils [35]. There are thus many opportunities for acquisition of NTM infection during ordinary daily activities, although the true incidence is not known. In contrast, nearly all transmission of Mycobacterium tuberculosis infection results from inhalation cough-generated aerosols from persons with active pulmonary TB. The annual risk of TB infection (ARTI) can be calculated from age-specific rates of tuberculin skin test reactivity; in most instances it is directly related to TB prevalence. Thus, although the ARTI may reach 4% in highly TB-endemic regions such as South Africa, it is as low as 0:01% in much of Northern Europe and North America [36][37][38]. These epidemiologic findings are consistent with the results of our mathematical model, and underscore the interplay of microbial and host biology in determining the relative contributions of reinfection and reactivation to mycobacterial pathogenesis.
Hence, while the screening for TB is necessary prior to initiating biologics, when considering NTM only a watchful monitoring during the treatment is recommended. This finding is particularly relevant, since it allows to avoid screening for NTM infection, which is complicated by the poor sensitivity of chest radiograph and more expensive and invasive techniques, such as chest computed tomography scan and/or bronchoscopy, should be used.
It is worth stressing that this methodology, being based on very standard stochastic procedures, has the advantage to hold beyond the test case of Mycobacteria. We hope that this test-case may shed light to future developments of this sideline approach in figuring out adverse events of biological therapies.   19). Here we used t~10years, consistently with clinical data. Notice that here parameters only count through the combination PN, and, again, N was thought as drug-independent. Moreover, in this case available data are extensive, so that the explicit number of overall treated patients N has been introduced. The average relative error on these parameters is &8%. doi:10.1371/journal.pone.0055017.t002