Mathematical model of oxygen, nutrient, and drug transport in tuberculosis granulomas

Physiological abnormalities in pulmonary granulomas–pathological hallmarks of tuberculosis (TB)–compromise the transport of oxygen, nutrients, and drugs. In prior studies, we demonstrated mathematically and experimentally that hypoxia and necrosis emerge in the granuloma microenvironment (GME) as a direct result of limited oxygen availability. Building on our initial model of avascular oxygen diffusion, here we explore additional aspects of oxygen transport, including the roles of granuloma vasculature, transcapillary transport, plasma dilution, and interstitial convection, followed by cellular metabolism. Approximate analytical solutions are provided for oxygen and glucose concentration, interstitial fluid velocity, interstitial fluid pressure, and the thickness of the convective zone. These predictions are in agreement with prior experimental results from rabbit TB granulomas and from rat carcinoma models, which share similar transport limitations. Additional drug delivery predictions for anti-TB-agents (rifampicin and clofazimine) strikingly match recent spatially-resolved experimental results from a mouse model of TB. Finally, an approach to improve molecular transport in granulomas by modulating interstitial hydraulic conductivity is tested in silico.


Introduction
Tuberculosis (TB) afflicts roughly one-third of the global population and causes approximately 1.5 million deaths annually [1].Its treatment often requires a lengthy drug therapy, in part because the dormant Mycobacterium tuberculosis bacilli hide within the core of hardened cellular masses in lungs called granulomas that offer a significant barrier to the transport of nutrients and therapeutic drugs.It is, thus, important to understand the flux of small molecules within granulomas, so that strategies might be devised to overcome transport limitations and hence improve treatment outcomes.
In order to reach the granuloma-contained bacilli, drugs must travel via the circulatory system to the lesion site, exit the granuloma vasculature (i.e., transvascular transport), and traverse the granuloma interior (i.e., interstitial transport) before encountering its target (i.e., cellular metabolism) [2].A related area of study wherein the path of circulating agents to and through abnormal masses has been thoroughly explored is the field of tumor transport [3,4].Indeed, TB granulomas and cancerous solid tumors are morphologically similar, resulting in similar transport limitations that have implications for therapeutic delivery and efficacy.As we have previously reported, a shared characteristic of these two cellular structures is an abnormal vasculature that is non-uniformly distributed within the tissue, resulting in avascular, diffusion-dominated regions, and thus, poor oxygen and drug delivery [5].
Within tumors, transvascular transport is compromised by the poorly formed and leaky tumor vasculature [6].This negatively impacts interstitial transport, exacerbated further by a lack of functional lymphatics, which in normal physiological settings serve a crucial role in homeostatic fluid balance within tissues [7].Impaired drainage by non-functional lymph vessels combined with abnormal, hyper-permeable blood vessels, results in an overall reduced pressure difference between the blood vessel and the interstitium [7], which is the driving force for transcapillary plasma exchange.This increases the interstitial fluid pressure (IFP), thereby reducing the effectiveness of the convective delivery of oxygen, nutrients, and drugs from the blood vessels within the core of the tumors, and hence providing a considerable barrier to effective delivery of anti-tumor agents.Although an absence of lymph vessels has not been confirmed within TB granuloma masses, we posit here that a similar IFP rise occurs within TB granulomas, as we have previously shown that granuloma blood vessels are structurally and functionally abnormal in the same manner as tumor blood vessels-and can even be "normalized" using the same targeted pharmacological approaches as for cancer [5].
The rise in IFP toward the core causes a radially outward convective flow of interstitial fluid velocity (IFV) oozing out at the periphery, thwarting the inward diffusion of drugs and nutrients [8,9].Interstitial diffusion, particularly of macromolecules, is further impeded by a dense cellular matrix.This has motivated many studies of tumor transport so that impediments to anti-cancer therapies can be understood and strategies developed to overcome them [10][11][12][13][14].Although IFP and IFV have not yet been directly measured within granulomas as they have been in tumors, we posit that a similar phenomenon occurs in TB, as we have found that other aspects including vasculature and hypoxia have been found to be analogous between the two lesion types in our prior experimental work [5].
The observable physiological abnormalities motivate a better understanding of the disease etiology and the underlying physicochemical properties that result in granuloma transport limitations.We first addressed this via a mathematical analysis of interstitial diffusion to quantitatively describe the emergence of hypoxic and necrotic granuloma regions as a direct result of limited oxygen availability [15].We now build upon these findings to include both convection and diffusion in the variably perfused regions in the granuloma microenvironment (GME).Because of the morphological similarities between tumors and granulomas, the substantial tumor transport modeling and experimental literature was leveraged to guide this approach; thus, results are shown for both abnormal masses.

Materials and methods
The description here is limited to the physiological basis (Fig 1) and main equations.The complete model derivation, assumptions, boundary conditions, parameter values, analytical solutions, and additional findings are provided in S1 Text.Results were plotted in Figs 2-7 for IFP (from Eq S30), IFV (from Eq S33), and chemical species concentrations (from Eq S14) from the main equations described below in Model Formulation and Assumptions and in S1 Text, using the built-in function "NDSolve" in Mathematica (Wolfram Research Inc., Champaign, IL).

Model formulation and assumptions
We assume here that the transcapillary exchange of oxygen is limited to an outer thoroughly vascularized shell [5,15].The rest of the GME is considered as completely avascular, thus allowing the assumptions of: 1) no extravasation, and 2) only diffusive (and no convective) transport therein (Fig 1).(In reality, the blood microvessel density declines more gradually [5], instead of an abrupt transition between well-perfused and inner regions assumed in our shell-core model; see S1 Fig) .The main equations governing the figures shown in the Results are as follows (see S1 Text for the governing equations and derivation).
Transport of oxygen and other small molecules in an idealized, spherical, vascularized granuloma can be written as a one-dimensional (radial) mass balance in the interstitial space (equivalent to Eq S11 in S1 Text) |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl }

Extravasation and Dilution ðEqIÞ
Where, for any chemical species j (e.g., oxygen, a nutrient, or a drug), its concentration is C j , which varies with respect to time (t) and radial location (r) within the idealized granuloma (Fig 1).Compared to our prior computational model (which only considered diffusion and reaction), we now include convection, extravasation, and dilution terms which we describe in detail here.On the left-hand side: 1) the accumulation term becomes 0 assuming steady-state (i.e., no accumulation of oxygen in the GME), and 2) the convection (second term on the left; v r is the radial interstitial fluid velocity) and diffusion (third term on the left; D j e is the effective diffusion coefficient of species j in the interstitial fluid) terms account for interstitial transport of species j.On the right-hand side: 3) the reaction term is assumed to follow Michaelis-Menten kinetics of cellular consumption (e.g., of oxygen, where k is the first-order rate constant for species consumption and K is the inverse of the half-saturation Michaelis-Menten constant) [15], while 4) the last term includes transvascular transport (i.e., extravasation) and plasma dilution (where L p is the hydraulic conductivity of the membrane-like blood vessel wall, a v is the volumetric vessel surface area in the tissue, p ve is the effective vessel pressure, p i is the interstitial fluid pressure, and C j,b is the bulk concentration of the species j in plasma).In dimensionless form (dimensionless terms defined in Eq S12 in S1 Text), the steady-state mass balance (in terms of dimensionless species concentration f over dimensionless radial distance ξ) reduces to (equivalent to Eq S14 in S1 Text) where Pe is the Pe ´clet number, i.e., the ratio of convective to diffusive transport, ν is the dimensionless radial interstitial velocity, ϕ 2 is the Thiele modulus, i.e., the ratio of reaction to diffusion rates, χ is the dimensionless Michaelis-Menten kinetic factor, α is the dimensionless modulus (i.e., granuloma size), ω is the dimensionless diffusion rate, and % is the relative extravasation/dilution rate.Explicit expressions for interstitial fluid pressure and velocity derive from the steady-state continuity equation, subject to suitable assumptions and boundary conditions (Eq S14 and S21) as described in the S1 Text.Upon testing multiple perfusion distribution scenarios (S1-S3 Figs), the idealized case of uniform perfusion was selected as the final form used here, given the striking similarity of its analytical solutions to the physiological reality of non-uniform perfusion (S3 Fig) .In dimensionless forms, the radial IFP and IFV are, respectively (equivalent to Eqs S30 and S33 in S1 Text) where ψ is dimensionless IFP, ν is dimensionless IFV, and α 0 is the dimensionless form of the outer granuloma radius (i.e., maximum dimensionless modulus).

Experimental data and parameter values
To the best of our knowledge, although IFP and IFV are well-studied in cancerous tumors, they have thus far not been theoretically investigated nor experimentally measured in TB granulomas.This means that any corroboration of the theoretical results must depend on corresponding experimental results for tumors.We predict TB granuloma pressure and velocity profiles in the Results, hypothesizing that granuloma properties are similar to those of tumors (Eq S39 and Eq S41 in S1 Text) based on our previous observations of morphological and functional similarities between cancerous tumors and TB granulomas in experiments [5].Using tumor parameter values, the dimensionless granuloma size (i.e., the modulus α) should be in the range of α 0 = 1−15 in our simulations, as a common granuloma diameter in the rabbit model is 2R 0 = 0.5−5 mm (although they can coalesce into larger lesions) [5,16].We compare these simulations to IFP/IFV data from tumor models [9].The oxygen consumption Michaelis-Menten parameters are utilized based on our previous work in rabbit TB granulomas [15].Drug delivery data from mouse models of TB [17] are compared to our simulations, with estimated tissue diffusion coefficients informed by anti-TB drug uptake in rats [18].All parameters are summarized in S1

Interstitial Fluid Pressure and Velocity Radial Profiles
We and others have experimentally and theoretically investigated the interstitial fluid pressure (IFP) in tumors [8,[19][20][21][22][23][24][25][26].We hypothesize here that there is an analogous IFP rise within TB granulomas, which IFP rises with distance into the granuloma, while effusive IFV follows an opposite trend and decreases towards the interior of the granuloma.We theorize in analogy to tumors [8] that i) the plasma that leaks out of the abnormal blood vessels as previously described [5] causes this rise in the IFP towards  the granuloma core and ii) this is exacerbated by a paucity of functional lymph vessels As IFP rises, the pressure drop declines until it eventually becomes zero.This is given by the Starling equation that accounts for the plasma flux across the microvessel walls between i) the effective vessel pressure (p v,e � p v −σ S (π v −π i )), where p v is the vascular pressure (mmHg), π v and π i are the vascular and interstitial osmotic pressures (mmHg), respectively, and σ S is the Staverman reflection coefficient) and ii) the IFP (p i , mmHg).This is especially apparent for larger granulomas (Fig 2A, α 0 >10), where the IFP approaches the effective vessel pressure rapidly with increasing distance towards the core.
As a result of the reducing pressure differential, the effusive (or filtration) flux of plasma, ϕ B , from the blood vessels (bearing oxygen, nutrients, and drugs) declines in the inward radial direction eventually becoming zero.Consequently, at this point, there is no further source for oxygen as well as other nutrients or small molecules/drugs, so that beyond this point the radial flux of these entities is entirely by their diffusion along with any interstitial (albeit compromised) fluid convection.Thus, IFV drops and stagnates in the core; this is especially apparent for large granulomas (Fig 3A, α 0 >10) where the IFV drops to zero rapidly with increasing distance into the granuloma.At the periphery of the granuloma (r !R 0 ) we may infer that, as in tumors [8], the IFV is directed outwards, which would further impede transport in the GME.

Convective zone thickness
As depicted schematically in Fig 1 , convection dominates in the vascularized, well-perfused outer rim of the GME, i.e., the region closest to the normal lung parenchyma.We can estimate the thickness of this rim (Δ = R 0 −R D , Fig 1 ) by assuming that within this shell, the pressure gradient becomes small enough that convective velocity approaches 0, representing the turning point at which diffusion dominates (i.e., this can be defined as the limited perfusion velocity; see S1 Text).Utilizing this assumption in Eq IV (equivalent to Eq S33 from the S1 Text), the dimensionless thickness of the convection zone, λ Δ , is obtained and plotted vs. granuloma size for varying limited perfusion velocities (Fig 4).For small granulomas (α 0 <5), the majority of the GME is predicted to be well-vascularized; this aligns with experimental observations [5,16].For 5<α 0 <10, the convective zone thickness drastically declines with increasing granuloma size until, at α 0 >10, it plateaus.This indicates that for large granulomas, the convection zone occupies less than half of the total mass thickness and the thickness of the convection zone becomes independent of total granuloma size.Utilizing our previously collected data for convective zone thickness from the rabbit TB model [15], in conjunction with theoretical parameter estimates from tumors (Eq S39 in S1 Text), we observe good agreement between our predicted and experimental results.

Oxygen and glucose concentration profiles
Expanding our previous modeling of avascular granulomas limited to diffusive transport combined with cellular consumption [15], our current model predicts oxygen concentration profiles while accounting for transcapillary exchange, plasma dilution, and interstitial convection and diffusion (Fig 5A).As expected, oxygen concentration drops quickly with distance into the granuloma core with increasing granuloma size.For α 0 >10, oxygen is predicted to be fully consumed rapidly within the GME; in reality, we know from prior measurements that necrosis emerges due to cell death at oxygen levels nearing-but not reaching-anoxia [5,16].Similar profiles result from the modeling of glucose transport (Fig 5B)-a slightly larger chemical species with poorer diffusivity (S1 Table ) than molecular oxygen-resulting in more rapid consumption of this nutrient particularly for α 0 >10, wherein the central GME is predicted to be devoid of glucose.

Overcoming interstitial transport barriers
The effect of tissue hydraulic conductivity on oxygen transport can be observed via in silico perturbation, whereby the simulated dimensionless oxygen concentration profile within the GME is shown for when the estimated base case parameter of tissue hydraulic conductivity K v (Eq S39 in S1 Text) is increased by a factor of 10 for small (α 0 = 3.5) and large (α 0 = 20) granulomas (Fig 6A).For both small and large granulomas, increasing hydraulic conductivity does not significantly affect oxygen concentration.However, in the large granuloma, the base value hydraulic conductivity value results in the glucose concentration dropping quickly to zero around a dimensionless radius value of 0.5 (Fig 6B).This implies that this nutrient is theoretically absent from the central GME under base tissue hydraulic conductivity.Significantly increasing tissue hydraulic conductivity (e.g., by drug treatment), results in glucose not being fully consumed before reaching the granuloma core.Indeed, a paucity of glucose is no longer predicted to emerge even in this very large simulated granuloma with the improvement in hydraulic conductivity, which suggests improved glucose distribution in the granuloma core.

Model comparisons to experimental drug delivery measurements in granulomas
Finally, this model can be applied not only to oxygen and nutrient delivery, but to anti-TB drug delivery as well.In a recent study [17], Kokesch-Himmelreich et al. measured drug delivery of rifampicin (RIF; molecular weight = 823 g/mol), clofazimine (CFZ; molecular weight = 473 g/mol), and other antibacterial agents using the spatially-resolved matrix-assisted laser desorption/ionization mass spectrometry imaging (MALDI-MSI) method in murine TB granulomas.We extracted and non-dimensionalized the spatially-resolved drug delivery data for RIF and CFZ in four granulomas (all R 0 <0.5 mm, α 0 <5), from this paper [17].Strikingly, our simple mass transport model accurately predicts the delivery of these two anti-TB agents, with low MSE values of ~0.01 for both drugs (Fig 7

Discussion
A detailed understanding of the drug transport barriers within granulomas can illuminate causes underlying the necessary prolonged treatment for TB, and suggest approaches to alleviate these physiological abnormalities of the GME that hinder transport.Thus, analogous modeling approaches from cancer research are applied here, based on the structural and morphological similarities between TB granulomas and tumors [5] despite the obvious differences between the two disease etiologies.Indeed, common modeling approaches for tumors that have been explored in TB granulomas by us and others include reaction-diffusion [15,27], continuum [28,29], agent-based [30][31][32], and multiscale [33][34][35] models, and vascular network/ angiogenesis models from tumors [36,37] are likely to be applied to TB granulomas in the future.The approximate analytical solution we obtained previously [15], for the case of diffusion-limited reaction, was able to predict the size of hypoxic and necrotic regions in good agreement with our experimental results from the rabbit TB model.Here, we extended our mathematical framework to include the effect of convective transport within the GME as a result of its abnormal vasculature, and to determine its impact on oxygen, nutrient, and drug delivery.

Model features, findings, and limitations
Transport phenomena including pressure rise, velocity profile, and interstitial diffusion, are dependent on the structural properties of the GME.We assume that the TB granuloma structure is similar to that of a tumor [21], which is comprised of three regions [5,15]: 1) the vascular region, limited to the periphery and devoid of functional lymphatic vessels; 2) the cellular region containing the cells; and 3) the interstitial (extra-cellular) space that contains an extracellular matrix network that imparts mechanical rigidity, bathed in interstitial fluid.We have shown previously that granuloma vessels are structurally and functionally abnormal, leading to inhomogeneous transport [5].In tumors, this aberrant vasculature compromises transvascular and interstitial convection due to high vessel permeability and an increase in IFP from plasma leakage and a lack of functional lymphatics.Together, these effects result in an outward convective interstitial velocity that opposes inward transport.By assuming that the GME recapitulates these transvascular and interstitial transport limitations seen in tumors based on our prior experimental evidence [5], the objective of our modeling effort was to predict these parameters and to indicate fruitful directions for future experimental investigations.
We provide here a comprehensive theoretical model of oxygen transport and reaction within a granuloma (or a tumor) that accounts for non-uniform vasculature, transcapillary exchange, plasma dilution, and interstitial convection and diffusion.Three limiting models of vasculature distribution were considered simultaneously for the first time: a model of uniform MVD distribution, a shell-core model, and a non-uniform MVD distribution model.Based on our analyses for IFP and IFV, we can conclude that the simpler case of uniform MVD distribution is adequate for predicting species transport and reaction.In contrast to earlier tumor IFP/ IFV models [8,19], our model accounts for plasma dilution (which is often excluded in the literature as negligible).Our mass balance formulation accounts for plasma extravasation from blood vessels not only as a source for oxygen, nutrients, and drugs, but as a diluent as well.Indeed, our model of interstitial transport was found to be in accord with tumor experiments.Predictions of convective zone thickness were found to be in agreement with experimental data from granulomas.The model was also utilized to theoretically investigate the effect of enhanced tissue hydraulic conductivity for overcoming transport barriers.Finally, our model accurately predicted the delivery of two anti-TB agents-rifampicin (first-line therapy) and clofazimine (second-line therapy)-with surprisingly good agreement given the simplicity of the mathematics describing the underlying transport phenomena.
In short, this model sheds light on the limitations of oxygen, nutrient, and drug transport within granulomas and tumors, and how such barriers might be overcome.For example, it predicts an IFP rise towards the granuloma core as a result of the known abnormal vasculature.This in turn results in a reduced or even stagnated IFV, compromising interstitial convective transport and thus depriving the granuloma core-where the bacilli hide-of oxygen, nutrients, and drugs.Our model predicts an absence of glucose in the granuloma core; this warrants experimental confirmation.The persistent bacteria in the hypoxic and necrotic core are able to utilize alternative carbon sources (e.g., lactate) [38]; thus, these metabolic alterations may present targetable vulnerabilities for therapeutic strategies.Modulating the interstitial hydraulic conductivity, e.g., with an agent targeting the interstitial components (such as losartan, a widely-prescribed, safe and inexpensive anti-hypertensive drug [39]), could prove to be an effective host-directed therapeutic strategy to modulate the GME and improve drug delivery and efficacy.
There are some limitations to the model.As described in detail in S1 Text, the non-uniform and/or shell-core models of perfusion better represent the physiological reality of blood vessel distribution in granulomas than the uniform perfusion model.However, the simpler limiting case of uniform perfusion advantageously allows for analytical solutions that 1) do not require artificially defining regions where convection or diffusion dominate, 2) are in good agreement when fitted with the non-uniform perfusion numerical solutions, and, most importantly, 3) accurately predict IFP/IFV experimental data.We also ignore cell membrane transport as a factor in the mass balance equation; however, we account for the cellular reaction of the species (i.e., oxygen consumption).Because we do not have experimental parameter measurements for TB granulomas (e.g., membrane and interstitial hydraulic conductivities) we apply tumor parameter values to the model given the similarities between granulomas and tumors, based on our previous experimental and computational observations regarding shared morphological and functional characteristics of these diseased masses (including vascular density, architecture, and perfusion) [5,15].We also do not consider variable oxygen/nutrient/drug uptake rates by different types of immune cells, e.g., macrophages vs. T cells, which we have recently demonstrated experimentally can contribute to heterogeneous drug distribution within the GME [40].Finally, because we do not consider the transport of large molecules here, we ignore any retardation factor.Indeed, because anti-TB agents often bind to proteins such as albumin [41], it may become necessary to model the unbound vs. bound drug fractions in future considerations.

Future directions
To the best of our knowledge, we provide here the first consideration of convective transport in the GME, and initial predictions for IFP and IFV profiles and their associated consequences in TB granulomas.To further support these findings, it should be experimentally confirmed whether-as in tumors-TB granulomas lack functional lymphatics.It should be noted that the first predictions of tumor IFP/IFV [8] preceded, yet accurately predicted, the first experimental measurements of these parameters by a number of years [42][43][44][45][46][47].We posit that future IFP/IFV measurements in TB granulomas will similarly confirm our predictions here.Continued comparison of tumors and TB granulomas via modeling may reveal novel similarities and differences between these two types of masses, particularly with regards to immune state and function.Furthermore, our findings support future testing of host-directed therapies that can modulate the GME to overcome transport barriers and improve treatment outcomes for this virulent and deadly disease.Indeed, following recent efforts by others to computationally optimize the implementation of multiple antibiotic treatments [48], future modeling efforts should provide rational basis for multi-drug dosing and scheduling.

Fig 1 .
Fig 1. Physiological basis for compartmentalized transport models in TB granulomas.Schematic depicting regions and consequences of compromised oxygen transport in idealized spherical granulomas, including 1) a vascularized region where convection dominates and plasma filtration from blood vessels occurs, and 2) an inner region lacking blood vessels where diffusion dominates, and hypoxia and necrosis result.(Adapted from [15]).https://doi.org/10.1371/journal.pcbi.1011847.g001

Fig 2 .
Fig 2. Granuloma IFP estimates and comparisons to tumor data.(A) Predicted dimensionless IFP profiles within granulomas from the uniform perfusion case for varying values of dimensionless granuloma size, α 0 = 1-15.(B) Fitting the theoretical IFP estimates (Eq III with a fitted modulus a 2 0 ¼ 24:4 [9], see S1 Text) to experimentally measured tumor IFP data (from human neuroblastoma tumor models grown in immunosuppressed rats, ~2 cm in diameter [9]) demonstrates the applicability of the uniform perfusion model to physiological IFP levels with a single fitted parameter.https://doi.org/10.1371/journal.pcbi.1011847.g002

Fig 6 .Fig 7 .
Fig 6.Effect of tissue hydraulic conductivity on oxygen and glucose delivery.(A) Oxygen and (B) glucose concentration profiles for base case parameter values (Eq S39, see S1 Text) of tissue hydraulic conductivity K v increased by a factor of 10 for small (α 0 = 3.5) and large (α 0 = 20) granulomas.https://doi.org/10.1371/journal.pcbi.1011847.g006 DiResta et al. measured IFP and IFV in human neuroblastoma tumor models in rats [9].Our predictive models of IFP (Fig 2B) and IFV (Fig 3B) align well with these experimental measurements (despite large error bars from the in-situ measurements) using a fitted value of α 0 = 4.95, analogous to S3 Fig. Thus, the simple analytical solution for the uniform perfusion case is strikingly accurate in predicting in vivo IFP/IFV levels, even in large lesions that are presumably non-uniformly perfused and likely contain hypoxic and necrotic regions.