A spatially heterogeneous network-based metapopulation software model applied to the simulation of a pulmonary tuberculosis infection

Tuberculosis (TB) is an ancient disease that, although curable, still accounts for over 1 million deaths worldwide. Shortening treatment time is an important area of research but is hampered by the lack of models that mimic the full range of human pathology. TB shows distinct localisations during different stages of infection, the reasons for which are poorly understood. Greater understanding of how heterogeneity within the human lung influences disease progression may hold the key to improving treatment efficiency and reducing treatment times. In this work, we present a novel in silico software model which uses a networked metapopulation incorporating both spatial heterogeneity and dissemination possibilities to simulate a TB infection over the whole lung and associated lymphatics. The entire population of bacteria and immune cells is split into a network of patches: members interact within patches and are able to move between them. Patches and edges of the lung network include their own environmental attributes which influence the dynamics of interactions between the members of the subpopulations of the patches and the translocation of members along edges. In this work, we detail the initial findings of a whole-organ model that incorporates distinct spatial heterogeneity features which are not present in standard differential equation approaches to tuberculosis modelling. We show that the inclusion of heterogeneity within the lung landscape when modelling TB disease progression has significant outcomes on the bacterial load present: a greater differential of oxygen, perfusion and ventilation between the apices and the basal regions of the lungs creates micro-environments at the apex that are more preferential for bacteria, due to increased oxygen availability and reduced immune activity, leading to a greater overall bacterial load present once latency is established. These findings suggest that further whole-organ modelling incorporating more sophisticated heterogeneities within the environment and complex lung topologies will provide more insight into the environments in which TB bacteria persist and thus help develop new treatments which are factored towards these environmental conditions. Electronic supplementary material The online version of this article (10.1007/s41109-018-0091-2) contains supplementary material, which is available to authorized users.

The Network package reflects the environment which is being simulated; this takes the form of a network of nodes and edges. The nodes of the network are Patches and each of these contain a subsection of the total system population. These sub-populations are divided into distinct compartments -each of which has a count reflecting how many members of that compartment exist within the patch at the current time. A member may only belong to one compartment at any given time, but may switch to a different compartment dependent on the system dynamics and the compartment paradigms chosen (i.e two compartments may model a real-world species in different states, such as infected or susceptible in epidemic dynamics). Each patch may also contain a set of attributes, which allow for the creation of heterogeneous environments, with different areas of the environment having different resources which influence the interactions and translocations that occur upon or between them. The patches are joined by instances of an Edge class. Each edge indicates some form of possible dispersal of members from one patch to another. Edges may also include attributes, indicating some difference in the environment that favours or hinders movement in a specific direction. Edges may be set as bi-directional or directed, implying free movement between two patches or movement in only one direction, respectively.

Events
The Events package is used to simulate the interactions between individuals within the patches and translocation of members between patches. Event modelling within the MetapopPy system uses the Gillespie Algorithm 1 to perform events asynchronously [1].
An event class is initialised with a reaction parameter (which can be interpreted as the rate of a single event occurrence in the specified unit time), and a state variable (which can be interpreted as the number of possible reactions in the network for the given event. The rate for the event is thus computed as the product of the reaction parameter and the state variable. When a change occurs on the network, each event must recalculate its rate by recalculating the state variable using the new population values. Each event class must define a function to calculate its state variable by interrogating the network: this function determines the number of possible reactions for the given event type based on the contents of the patches and edges of the network. Each event instance is initialised with a reaction parameter (which can be interpreted as the rate of a single event occurrence in the specified unit time) and the event calculates its rate by multiplying the reaction parameter and the state variable.
Each event class must also define a perform function, which determines what changes are made to the subpopulations of the network when the event occurs. These functions are used within the Dynamics package. An event is classed at occurring at only one patch, but may update multiple patches. For example, movement from patch A to patch B is classed as occurring at patch A and depends only on the values of the population there (and the edges attached to patch A), but when performed, the event updates the population of both patches.

Dynamics
The Dynamics package joins a network with a series of events that enact upon it. When a simulation of the model is run, the network is seeded with initial conditions: a set of subpopulation compartment counts (determined by the user) and attributes for all patches. The simulation time is set initially at zero. MetapopPy then uses an asynchronous time modelling system following the Gillespie algorithm: each event calculates its rate as specified above and an event is chosen probabilistically based upon these rates, with a greater rate indicating a greater likelihood of being chosen. A timestep, τ , for this event to occur is then chosen using Equation 1 (as per [1]), with a being the total sum of the rate of all events, and r being a randomly generated number in the range (0,1).
The chosen event is then performed, and updates the network in its defined manner. The simulation time is then incremented by τ and the process repeats, until a set time-limit is exceeded or there is no possibility of any event occurring (i.e. a=0). The sub-population values of all patches are recorded at a userdefined time interval throughout. When a patch is updated by an event, other events need to recalculate their state variables and their rates. To recalculate the state variable contribution of all patches to all events would be inefficient, as some values may not have changed. Therefore, an UpdateHandler class is included to propagate updates efficiently. This class tracks which events occur at which patch, and also which compartments are used to calculate the state variable contribution for each event. When the compartments of a patch are updated, the update handler is triggered and only events which a) occur at the patch, and b) are dependent on the amended compartments, recalculate their state variable contribution from the amended patch. By tracking these dependencies and utilising them in this manner, the computation needed after each patch update is drastically reduced.

TBMetapopPy Events
TBMetapopPy contains a number of events to simulate the growth of bacteria and the interactions between the immune cells of the body and the bacteria. These events are described in this section, and listed in tables 1, 2, 3 and 4, along with their respective reaction parameters, state variable equations and outcomes.

Bacteria events
An initial population of bacteria is deposited in the lung. M I members may become overwhelmed by their internal population of B IM bacteria, causing them to 'burst' (event 34), destroying the M I member and releasing the average internal number of bacteria back out. We assume that the internal conditions of the macrophage are a stressful environment, and therefore all bacteria are released into the B ED compartment, having become dormant to survive inside the immune cell.

Dendritic cell events
Dendritic cells act similarly to macrophages, and are recruited (event 18) and die (event 16) at constant rates, with recruitment of DC I members being increased by the presence of extracellular bacteria (event 19).
Contact between DC I and extracellular bacteria may cause uptake of a bacterium (events 13 and 14), converting DC I to DC M and the bacterium to B ID (we assume dendritic cells do not destroy bacteria). DC M members may translocate (event 15) to the LymphPatches along LymphEdges for antigen presentation.

T-cell events
A population of T N members exists in the LymphPatch instances at the start of simulation. Like other immune cells, these cells die (event 42) and are replaced by newly recruited cells (event 36). As the infection spreads and APCs (M I and DC M ) transfer to the lymphatics, these cells may trigger the naïve population to activate (events 39 and 40), converting T N into T A . These T A members may then migrate back to the lung (event 41) along a BloodEdge, with the choice of edge linked to the perfusion value, resulting in easier travel to areas with greater perfusion.
Activated T-cells have two main functions in the model: the first is to trigger activation of macrophages (converting M R to M A , event 29). Activated macrophages will naturally revert to a resting state (event 30) if the t-cell presence is too low.
Secondly, T-cells cause the cytotoxic destruction of infected macrophages (event 35). The M I macrophage is destroyed, and likewise any B IM bacteria inside are also destroyed. These functions may be performed at both classes of patch.