Computational Modeling of Microabscess Formation

Abstract

Bacterial infections can be of two types: acute or chronic. The chronic bacterial infections are characterized by being a large bacterial infection and/or an infection where the bacteria grows rapidly. In these cases, the immune response is not capable of completely eliminating the infection which may lead to the formation of a pattern known as microabscess (or abscess). The microabscess is characterized by an area comprising fluids, bacteria, immune cells (mainly neutrophils), and many types of dead cells. This distinct pattern of formation can only be numerically reproduced and studied by models that capture the spatiotemporal dynamics of the human immune system (HIS). In this context, our work aims to develop and implement an initial computational model to study the process of microabscess formation during a bacterial infection.

1. Introduction

The immune system is one of the most important and complex system of our organism. Despite great advances in recent years that shed light on its understanding and unravel the underlying key mechanisms behind its functions, there are still many functions of the human immune system (HIS) that are not well understood. Computational models of HIS dynamics can contribute to a better understanding of the relationship between cells and molecules of the HIS.

In this study, we developed a mathematical model of some cells and molecules of the HIS to reproduce the spatiotemporal dynamics of the initial formation of microabscesses during an immune response to a bacteria.

To reproduce these dynamics, we introduce a mathematical model composed of a system of partial differential equations (PDEs) that extends our previous models [1, 2] and defines the dynamics of representative cells and molecules of the HIS during the immune response to a bacteria. The model presented is descriptive, mechanistic, and deterministic; therefore, it enables the understanding of how different complex phenomena, structures, and elements interact during an immune response. In addition, the model’s parameters reflect the physiological features of the system, making the model appropriate for general use.

The remainder of this paper is organized as follows. First, the necessary biological background is presented. Next, related works are briefly discussed. This exposition is followed by a description of the mathematical model proposed in this work and the numerical scheme used to implement it. Then simulation results obtained from the proposed model are discussed, and, finally, our conclusions and plans for future work are presented.

2. Biological Background

The initial response of the host to a diverse array of biological stressors including bacterial infection, burns, trauma, and invasive surgery is an inflammatory response. Despite the growing understanding of the cellular and molecular mechanisms of inflammation, the complexity of the inflammatory response has challenged therapeutic development [3, 4]. A key reason for this conundrum has been speculated to be the difficulty of predicting the impact of manipulating individual components of the highly complex, nonlinear, and redundant inflammatory response [5]. Thus, progress would require a greater understanding of how components are organized. This makes systems biology based approaches appealing [6].

Most inflammatory reactions begin as a suppurative or purulent exudation process when the first line of cellular defense, the neutrophils, accumulate in the area. A suppurative process is characterized by the presence of pus (neutrophils mixed with cellular debris). Classically, there are three requisites of suppuration:(i)presence of neutrophils that release proteolytic enzymes,(ii)necrosis of some types,(iii)liquefaction.

This suppurative process may lead to the formation of microabscesses. A microabscess is a localized collection of dead cells, body fluids, microbes, and other cells of the HIS. The process of formation of a microabscess begins when a cell of the HIS encounters bacteria and warn other cells that there is a stranger in the host. Its “warn” is in the form of a class of biochemicals called cytokines, which beckon other HIS cells to come to the point of infection and surround the enemy.

Most often all this goes unnoticed because the first few immune system cells phagocytize (engulf and digest) the invaders and the battle is finished. But every so often, an invader has a trick to escape the immune response and cannot be killed by the phagocytes. Those bacteria continue to grow and to spew out whatever they do. More and more immune system cells, mainly neutrophils, congregate at the infection site trapping the pathogens in the center. If this microabscess is close enough to the body surface, it can be seen as a blob of pus under the skin. When a microabscess is well developed, it has a wall or capsule of fibrous connective tissue separating it from the surrounding tissue, helping to prevent any microbes present in the microabscess from spreading to other areas of the body. Thus, microabscesses can be considered as a natural strategy used to fight against infection.

Microabscesses are found in many different diseases, for instance, the authers in [7–10] present animal studies detailing the formation of liver microabscess and microabscess by different types of infections. Epidermal microabscess formation by neutrophils was also evaluated in [11–14]. Infection of the heart by bacteria (bacterial myocarditis [15]) or by viruses (viral myocarditis [16]) is also correlated with microabscess formation by neutrophils. The interaction between tumor cells and inflammatory cells plays an important role in cancer initiation and progression and was investigated in [17] for the case of tumor-infiltrating neutrophils in pancreatic neoplasia, where the pattern of microabscess formation by neutrophils was reported once again.

3. Related Work

This section presents and discusses other mathematical and computational models of the immune response. Essentially, two distinct approaches are used: agent-based models and ordinary differential equations (ODEs). These models have some features in common with our model. All models include representative cells and molecules of the innate immune system. For example, neutrophils, macrophages, and proinflammatory cytokines are modeled in the majority of models. Some models as our model consider the important interactions between endothelial cells, tissue cells, and cytokines. Despite some similarities with our model, none of the works focus on modeling microabscesses.

4. Models Based on Agents

In [18, 19], it was developed an agent-based model of the dynamics of some cells, such as polymorphonuclear leukocytes (PMNs) and mononuclear cells and molecules, such as TNF-a and IL-1, during the initial inflammatory response in the interface endothelium/blood at the capillary level. Some characteristics of the model are as follows:(i)all the cells are represented as agents whose behavior is close to the real;(ii)it considers the interactions between endothelial cells and circulating inflammatory cells at the blood/blood vessel-lining interface;(iii)the initial injury number (IIN) defines the number of tissue cells that are dead initially;(iv)the oxygen concentration is one important variable in the model;(v)the total tissue damage is represented as a deficit in the oxygen variable;(vi)the injury state of an endothelial cell depends on the available oxygen concentration; (vii)proinflammatory mediators and endothelial cells surface adhesion molecules are modeled by state variables;(viii)it considers a generic pathogen that causes the infection.

The work aims to reproduce the time course of the early inflammatory response associated with the Systemic Inflammatory Syndrome Response (SIRS)/Multiple Organ Failure (MOF) from massive trauma or large exposure to endotoxin. The objective of the simulations were to compare the results with the soluble TNF-receptor experiment [20] where the soluble TNF-receptor is tested as a therapeutic treatment for the sepsis. The author claims that his results generally replicate the results of several large-scale clinical trials of cytokine-directed antimediator agents.

In [6], an agent-based modeling (ABM) framework is proposed to study the nonlinear dynamics of acute inflammatory responses to LPS. Their work uses an agent-based approach to elucidate molecular interactions involved in the NF-𝜅𝛽 signaling pathway, coupled with the spatial orientation of various inflammation specific molecules and cell populations such as macrophages and T-helper cells. In their model, the propagation of LPS signaling across the system is studied by considering the coupling between extracellular signals and transcriptional response. Some relevant characteristics considered in the model are as follows:(i)molecular interactions, (ii)cellular heterogeneity,(iii)LPS/TLR4 signal transduction pathway,(iv)transcriptional response.

The proposed in silico model is evaluated through its ability to successfully reproduce a self-limited inflammatory response as well as a series of scenarios: a persistent (non)infectious response or innate immune tolerance and potentiation effects followed by perturbations in intracellular signaling molecules and cascades.

5. Models Based on ODEs

The model of [21] studies immunomodulatory strategies for treating cases of severe sepsis. They introduced and evaluated the concept of conducting a randomized clinical trial in silico based on simulated patients generated from a mechanistic mathematical model of bacterial infection, acute inflammatory response, global tissue dysfunction, and a therapeutic intervention. Trial populations are constructed to reflect heterogeneity in bacterial load and virulence as well as propensity to mount and modulate an inflammatory response. They constructed a cohort of 1000 trial patients submitted to therapy with one of three different doses of a neutralizing antibody directed against tumor necrosis factor (anti-TNF) for 6, 24, or 48hrs. Their focus was to assess the feasibility of using differential equation models to improve the design of clinical trials. This paper replicates in silico the general findings from actual clinical trials—that it is very difficult to design a treatment strategy that is effective over a broad range of sepsis patients.

6. Hybrid Models

In [22], a hybrid model that coupled an agent-based model [18, 19] and a system dynamics/differential equation model [21] was created using the System Dynamics tool within Netlogo [23]. They developed an interface between the agent-based and system dynamics models. The area of initial infection and subsequent interactions at this point were simulated with the agent-based model and certain aspects primarily related to the production and life cycle of circulating inflammatory cells was modeled with the system dynamics model. They implemented a simple system dynamics model that focused on systemic polymorphoneutrophils (PMN) production, maturation, sequestration, and release. Some of the flows inside the system dynamics model were influenced by the conditions within the agent-based model. The primary role of the system dynamics model was to manifest a delay between the elevation of the cytokines in the tissue and the increase in PMNs in the circulating blood. The objective of the work was to reproduce some of the results of the agent-based model [18, 19] using the hybrid model.

7. Mathematical Model

Our main objective is to develop a parameterized mathematical model of the human innate immune system that simulates the immune response occurring in a generic tissue. To achieve this goal, we first build a model of the immune response to LPS [1, 2]. In this work, we extend this model to reproduce the spatiotemporal dynamics of a bacterial infection and the process of microabscess formation.

The mathematical model simulates the temporal and spatial behavior of bacteria (𝐵), dead bacteria (BD), macrophages, neutrophils (𝑁), apoptotic neutrophils (ND), proinflammatory cytokines (CH), healthy tissue cells (HT), and dead tissue cells (TD). Macrophages are present in two states of readiness: resting (RM) and hyperactivated (AM). We must stress that the equations modeling proinflammatory cytokines are generic in the sense that they model the role of distinct cytokines taking part in the inflammatory process. Equation parameters can be adjusted to model the role of a specific proinflammatory cytokine.

The relationships among all of the model’s components are described next. Neutrophils, resting macrophages, and active macrophages phagocytose the bacteria. The neutrophils then undergo apoptosis, which may or may not be induced by the phagocytosis process. In this different state, apoptotic neutrophils cannot perform phagocytosis or produce proinflammatory cytokines; as a result, apoptotic neutrophils are eliminated from the body after being phagocytosed by active macrophages. Apoptotic neutrophils will die after a period of time, releasing cytotoxic granules and degradation enzymes in the medium that cause tissue damage destroying healthy tissue cells. Active neutrophils and bacteria also cause tissue damage by producing toxic products that are not modeled here. The infection site is “cleaned” by resting and active macrophages that do the phagocytosis of dead tissue cells. Healthy tissue cells in contact with bacteria, neutrophils, and active macrophages produce proinflammatory cytokines. The proinflammatory cytokines increase the permeability of the blood vessels; consequently, more neutrophils and monocytes are recruited to the infected tissue. In addition, the proinflammatory cytokines act as a chemoattractant substance to the resting macrophages, active macrophages, and neutrophils.

Below, we provide the equations derived from the model. Equation (1) provides the bacteria differential equation: 𝜕𝐵𝜕𝑡=𝑟_𝐵⋅𝐵⋅𝑔(𝑤)−𝜇_𝐵⋅𝐵−𝜆_𝑁|𝐵⋅𝑁⋅𝐵−𝜆_RM|𝐵⋅RM⋅𝐵−𝜆_AM|𝐵⋅AM⋅𝐵+𝐷_𝐵(),𝐵⋅dif𝐵,𝑤𝑥,𝑦,0=𝐵₀,𝜕𝐵(⋅,𝑡)𝜕𝑛_𝜕Ω=0. (1)

In this equation, 𝑟_𝐵⋅𝐵⋅𝑔(𝑤) denotes the reproduction term of the bacteria, where 𝑟_𝐵 is the rate of reproduction and 𝑔(𝑤) is a function of the total density of cells 𝑤 in a discretized area of the two-dimensional space at a specific time step. The 𝑤 variable is defined as 𝑤;𝑥,𝑦,𝑡=𝐵𝑥,𝑦,𝑡+BD𝑥,𝑦,𝑡+𝑁𝑥,𝑦,𝑡+ND𝑥,𝑦,𝑡+RM𝑥,𝑦,𝑡+AM𝑥,𝑦,𝑡+CH𝑥,𝑦,𝑡+HT𝑥,𝑦,𝑡+TD𝑥,𝑦,𝑡 (2)𝜇_𝐵⋅𝐵 denotes the decay of bacteria, where 𝜇_𝐵 is the rate of decay. 𝜆_𝑁|𝐵⋅𝑁⋅𝐵 denotes the phagocytosis of bacteria by neutrophils, where 𝜆_𝑁|𝐵 is the rate of this phagocytosis. 𝜆_RM|𝐵⋅RM⋅𝐵 denotes the phagocytosis of bacteria by resting macrophages, where 𝜆_RM|𝐵 is the rate of this phagocytosis. 𝜆_AM|𝐵⋅AM⋅𝐵 denotes the phagocytosis of bacteria by active macrophages, where 𝜆_AM|𝐵 is the rate of this phagocytosis. 𝐷_𝐵⋅dif(𝐵,𝑤) denotes bacteria diffusion, where 𝐷_𝐵 represents the diffusion coefficient and dif(𝐵,𝑤) is calculated in the following way: dif(𝐵,𝑤)=∇⋅𝑔(𝑤)∇𝑓(𝑤)𝐵−𝑓(𝑤)𝐵∇𝑔(𝑤). (3) The 𝑓 function models the probability of a cell being pushed from a site due to the pressure exerted by neighboring cells [24, 25]. This population pressure is modelled by a Hill equation [26]. It increases with the total density of cells 𝑤 occupying the same position in space and has a saturation in a high density of cells.

The 𝑓 function is defined as 𝑓(𝑤)𝑤=1+𝛼.𝛽+𝑤 (4)𝛼 and 𝛽 are constant values.

The 𝑔 function returns the percentage of free space in a discretized area of the two-dimensional space and its use is motivated by some important biological concepts such as quorum sensing/volume sensing [24, 25, 27–29]. The idea is that cells have a set of cell density sensing mechanisms and changes its behavior in crowded regions. In the context of our model it is used to limit the density of cells that occupy a discretized area of our two-dimensional domain. The 𝑔 function is defined as 𝑔(𝑤)𝑤=1−.total (5) The variable total represents the maximum density of cells that fits in a discretized area of the tissue.

The differential equation corresponding to dead bacteria (BD) is given as follows:𝜕BD𝜕𝑡=𝜇_𝐵⋅𝐵+𝜆_𝑁|𝐵⋅𝑁⋅𝐵+𝜆_RM|𝐵⋅RM⋅𝐵+𝜆_AM|𝐵⋅AM⋅𝐵−𝜆_AM|BD⋅AM⋅BD−𝜆_RM|BD⋅RM⋅BD+𝐷_BD⋅dif(BD,𝑤),BD𝑥,𝑦,0=𝐵₀,𝜕BD(⋅,𝑡)𝜕𝑛_𝜕Ω=0. (6) Here, note that 𝜇_𝐵⋅𝐵, 𝜆_𝑁|𝐵⋅𝑁⋅𝐵, 𝜆_RM|𝐵⋅RM⋅𝐵 and 𝜆_AM|𝐵⋅AM⋅𝐵 were defined previously. 𝜆_AM|BD⋅AM⋅BD denotes the phagocytosis of dead bacteria by active macrophages, where 𝜆_AM|BD is the rate of phagocytosis. 𝜆_RM|BD⋅RM⋅BD denotes the phagocytosis of dead bacteria by resting macrophages, where 𝜆_RM|BD is the rate of phagocytosis. 𝐷_BD⋅dif(BD,𝑤) denotes dead bacteria diffusion, where 𝐷_BD represents the diffusion coefficient and the function dif was defined previously.

The differential equation corresponding to the resting macrophage (RM) is given as follows:𝑃_RM=𝑃^max_RM−𝑃^min_RM⋅CHCH+keq_CH+𝑃^min_RM,source_RM=𝑃_RM⋅𝑀^max−(RM+AM),𝜕RM𝜕𝑡=−𝜇_RM⋅RM−𝜆_𝐵|RM⋅𝐵⋅RM+source_RM⋅𝑔(𝑤)+𝐷_RM⋅dif(RM,𝑤)−𝜒_RM⋅chemotaxis(RM,CH,𝑤),RM𝑥,𝑦,0=RM₀,𝜕RM(⋅,𝑡)𝜕𝑛_𝜕Ω=0. (7)𝑃_RM denote the increase in endothelium permeability and its effects on monocyte extravasation. The permeability of blood vessel endothelium is modeled by a Hill equation [26], which also has been used to model drug dose-response relationships [30]. The idea is to model the increase in the permeability of the endothelium in accordance with the number of proinflammatory cytokines deposited on the endothelium.

The calculation of 𝑃_RM involves the following parameters: (a) 𝑃^max_RM, the maximum endothelium permeability induced by the proinflammatory cytokine; (b) 𝑃^min_RM, the minimum endothelium permeability induced by the proinflammatory cytokine; (c) keq_CH, the number of proinflammatory cytokines that exert 50% of the maximum effect on permeability. source_RM⋅𝑔(𝑤) denotes the source term of macrophages, which is related to the number of monocytes that will enter into the tissue from the blood vessels. This number depends on the endothelium permeability 𝑃_RM and on the number of monocytes appearing in the blood (𝑀^max).

𝜇_RMRM denotes resting macrophage apoptosis, where 𝜇_RM is the apoptosis rate. 𝜆_𝐵|RM⋅𝐵⋅RM denotes the activation of resting macrophages, where 𝜆_𝐵|RM is the rate of activation. 𝐷_RM⋅dif(RM,𝑤) denotes resting macrophage diffusion, where 𝐷_RM represents the diffusion coefficient and the function dif was defined previously. 𝜒_RM⋅chemotaxis(RM,CH,𝑤) denotes resting macrophage chemotaxis, where 𝜒_RM is the chemotaxis rate and chemotaxis(RM,CH,𝑤) is calculated in the following way: chemotaxis(RM,CH,𝑤)=∇⋅RM𝑔(𝑤)𝑓(𝑤)∇CH. (8) The differential equation corresponding to the active macrophage (AM) is given as follows:𝜕AM𝜕𝑡=−𝜇_AM⋅AM+𝜆_𝐵|RM⋅𝐵⋅RM+𝐷_AM⋅dif(AM,𝑤)−𝜒_AM⋅chemotaxis(AM,CH,𝑤),AM𝑥,𝑦,0=AM₀,𝜕AM(⋅,𝑡)𝜕𝑛_𝜕Ω=0. (9) Here, note that 𝜆_𝐵|RM⋅𝐵⋅RM was defined previously. Above, 𝜇_AM⋅AM, 𝐷_AM⋅dif(AM,𝑤), and 𝜒_AM⋅chemotaxis(AM,CH,𝑤) denote the active macrophage apoptosis, diffusion, and chemotaxis, respectively, whereas 𝜇_AM, 𝐷_AM, and 𝜒_AM are the apoptosis rate, diffusion coefficient, and chemotaxis rate, respectively.

The neutrophil differential equation (𝑁) is given as follows:𝑃_𝑁=𝑃_𝑁^max−𝑃_𝑁^min⋅CHCH+keq_CH+𝑃_𝑁^min,source_𝑁=𝑃_𝑁⋅𝑁^max−𝑁,𝜕𝑁𝜕𝑡=−𝜇_𝑁⋅𝑁−𝜆_𝐵|𝑁⋅𝐵⋅𝑁+source_𝑁⋅𝑔(𝑤)+𝐷_𝑁⋅dif(𝑁,𝑤)−𝜒_𝑁(),𝑁⋅chemotaxis𝑁,CH,𝑤𝑥,𝑦,0=𝑁₀,𝜕𝑁(⋅,𝑡)𝜕𝑛_𝜕Ω=0. (10) In this equation, 𝑃_𝑁 denotes the increase in endothelium permeability and its effects on neutrophil extravasation. In the top equation, 𝑃_𝑁^max is the maximum endothelium permeability induced by proinflammatory cytokines, 𝑃_𝑁^min is the minimum endothelium permeability induced by proinflammatory cytokines, and keq_CH is the number of proinflammatory cytokines that exert 50% of the maximum effect on endothelium permeability.

Here, 𝜇_𝑁⋅𝑁 denotes neutrophil apoptosis, where 𝜇_𝑁 is the rate of apoptosis. 𝜆_𝐵|𝑁⋅𝐵⋅𝑁 denotes the neutrophil apoptosis induced by phagocytosis, where 𝜆_𝐵|𝑁 represents the rate of this induced apoptosis. source_𝑁⋅𝑔(𝑤) represents the source term of neutrophil, that is, the number of neutrophils entering the tissue from the blood vessels. This number depends on the endothelium permeability (𝑃_𝑁) and on the number of neutrophils in the blood (𝑁^max).

𝐷_𝑁⋅dif(𝑁,𝑤) denotes neutrophil diffusion, where 𝐷_𝑁 represents the diffusion coefficient and the function dif was defined previously. 𝜒_𝑁⋅chemotaxis(𝑁,CH,𝑤) denotes neutrophil chemotaxis, where 𝜒_𝑁 is the chemotaxis rate and chemotaxis(𝑁,CH,𝑤) was defined previously.

The differential equation corresponding to the apoptotic neutrophil (ND) is given as follows:𝜕ND𝜕𝑡=𝜇_𝑁⋅𝑁+𝜆_𝐵|𝑁⋅𝐵⋅𝑁−𝜆_ND|AM⋅ND⋅AM−𝜇_NDND+𝐷_ND⋅dif(ND,𝑤),ND𝑥,𝑦,0=ND₀,𝜕ND(⋅,𝑡)𝜕𝑛_𝜕Ω=0. (11) Here, note that 𝜇_𝑁⋅𝑁 and 𝜆_𝐵|𝑁⋅𝐵⋅𝑁 were defined previously, whereas 𝜆_ND|AM⋅ND⋅AM denotes the phagocytosis of the apoptotic neutrophil carried out by active macrophages, and 𝜆_ND|AM is the rate of this phagocytosis. 𝜇_NDND denotes the neutrophil necrosis, where 𝜇_ND is the rate of necrosis. 𝐷_ND⋅dif(ND,𝑤) denotes apoptotic neutrophil diffusion, where 𝐷_ND represents the diffusion coefficient and the function dif was defined previously.

The differential equation for the proinflammatory cytokine (CH) is given in as follows:𝜕CH𝜕𝑡=−𝜇_CH⋅CH+𝛽_CH|𝑁⋅𝑁⋅𝐵+𝛽_CH|AM⋅AM⋅𝐵+𝛽_CH|HT⋅HT⋅𝐵⋅𝑔(𝑤)+𝐷_CH⋅dif(CH,𝑤)CH𝑥,𝑦,0=CH₀,𝜕CH(⋅,𝑡)𝜕𝑛_𝜕Ω=0. (12) In this equation, 𝜇_CHCH denotes the proinflammatory cytokine decay, where 𝜇_CH is the decay rate. 𝛽_CH|𝑁⋅𝑁⋅𝐵 denotes the proinflammatory cytokine production by the neutrophils, where 𝛽_CH|𝑁 is the production rate. 𝛽_CH|AM⋅AM⋅𝐵 denotes the proinflammatory cytokine production by active macrophages, where 𝛽_CH|AM is the production rate. 𝛽_CH|HT⋅HT⋅𝐵 denotes the proinflammatory cytokine production by healthy tissue cells in contact with bacteria, where 𝛽_CH|HT is the production rate. 𝐷_CH⋅dif(CH,𝑤) denotes the proinflammatory cytokine diffusion, where 𝐷_CH represents the diffusion coefficient and the function dif was defined previously.

The differential equation corresponding to the healthy tissue (HT) is given as follows:𝜕HT𝜕𝑡=−𝜇_NDND−𝜆_𝐵|HT⋅𝐵⋅HT,HT𝑥,𝑦,0=HT₀,𝜕HT(⋅,𝑡)𝜕𝑛_𝜕Ω=0. (13)𝜇_NDND denotes the tissue damage caused by the release of toxic products from necrotic neutrophils. 𝜆_𝐵|HT⋅𝐵⋅HT denotes the tissue damage caused by bacteria, where 𝜆_𝐵|HT is the rate of damage.

The differential equation corresponding to the dead tissue (TD) is given as follows:𝜕TD𝜕𝑡=𝜇_NDND+𝜆_𝐵|HT⋅𝐵⋅HT−𝜆_RM|TD⋅RM⋅TD−𝜆_AM|TD⋅AM⋅TD,TD𝑥,𝑦,0=TD₀,𝜕TD(⋅,𝑡)𝜕𝑛_𝜕Ω=0. (14)𝜇_NDND and 𝜆_𝐵|HT⋅𝐵⋅HT were defined previously. 𝜆_RM|TD⋅RM⋅TD denotes the phagocytosis of dead tissue cells by resting macrophages, where 𝜆_RM|TD is the rate of phagocytosis. 𝜆_AM|TD⋅AM⋅TD denotes the phagocytosis of dead tissue cells by active macrophages, where 𝜆_AM|TD is the rate of phagocytosis.

The mathematical model presented here introduced some modifications to our previous model [2] with the aim to reproduce the microabscess formation. We included equations for the dynamics of the tissue to take into account some effects of infection such as tissue damage and production of cytokines by tissue cells. We also replaced the LPS equation in our previous model [2] by the bacteria equation with a term for reproduction of bacteria. Besides we modified the calculus of the diffusion and chemotaxis terms [24] (a) to limit the number of cells that are allowed to stay at the same time in the same area of the domain and (b) to reduce the efficiency of the diffusion and chemotaxis processes in overcrowded regions. More specifically, the method we implemented incorporates the following general mechanisms which may lead to dispersal of the population [24].(i)Population pressure: we assume that a high cell density results in increased probability of a cell being pushed from a site, for example, due to the pressure exerted by neighboring cells. This is achieved phenomenologically with the 𝑓 function of our mathematical model and the changes in diffusion and chemotaxis calculus. (ii)Limited space: here we assume that no more cells can enter a site above a total cell density. In our model, this is achieved with the incorporation of the 𝑔 function in the diffusion and chemotaxis calculations. (iii)Gradient detection: cells may detect and respond to a local gradient in the cell density and as a consequence cells can move to higher concentrations of the attractant substance.

8. Implementation

The numerical method used to solve the mathematical model was the Finite Difference Method [31], a method commonly used in the numeric discretization of PDEs.

A complex part of the resolution of the PDEs is the resolution of the convective term, the chemotaxis term. The development of numerical methods to approximate convective terms (in most cases not linear) have been subject of intense researches [32–35].

Different numerical approaches have been proposed for the discretization of the chemotaxis term [36, 37]. Our implementation is based on the finite difference method for the spatial discretization and the explicit Euler method for the time evolution. The discretization of the chemotaxis term (∇⋅(𝜒_𝑁𝑁∇CH)) uses the First-Order Upwind scheme [38]. Therefore, the precision of our numerical implementation is first-order in time (explicit Euler) and first-order in space (upwind scheme). The upwind scheme discretizes the hyperbolic PDEs through the use of differences with bias in the direction given by the signal of the characteristics’ speeds. The upwind scheme uses an adaptive or solution-sensitive stencil to numerically simulate more precisely the direction of information propagation.

In two-dimension, the upwind scheme approximates the chemotaxis term as the sum of the flux in the 𝑥 direction (res𝑋) with the flux in the 𝑦 direction (res𝑌). res𝑋 is the sum of the flux_left at the point 𝑥−delta𝑋/2 with the flux_right at the point 𝑥+delta𝑋/2 and res𝑌 is the sum of the flux_up at the point 𝑦−delta𝑌/2 with the flux_down at the point 𝑦+delta𝑌/2 in Algorithm 1.

float chemotaxis (float** 𝑛, float **ch, int 𝑥, int 𝑦, float **𝑤) {  flux_left = 0;  flux_right = 0;  if (𝑥 > 0)  {     if ((ch[𝑦][𝑥]−ch[𝑦][𝑥−1])>0)     {      flux_left = − (ch[𝑦][𝑥] – ch[𝑦][𝑥−1])*𝑛 [𝑦][𝑥−1] *𝑔(𝑤[𝑦][𝑥−1])*𝑓 (𝑤[𝑦][𝑥−1])/delta𝑋;    }    else    {      flux_left = − (ch[𝑦][𝑥] – ch[𝑦][𝑥−1])*𝑛[𝑦][𝑥] *𝑔(𝑤[𝑦][𝑥])*𝑓(𝑤[𝑦][𝑥])/delta𝑋;    }  }  if (𝑥 < (size − 1))  {    if ((ch[𝑦][𝑥+1] – ch[𝑦][𝑥]) > 0)    {      flux_right = (ch[𝑦] [𝑥+1] – ch[𝑦][𝑥])*𝑛[𝑦] [𝑥]*𝑔(𝑤[𝑦][𝑥])*𝑓(𝑤[𝑦][𝑥])/delta𝑋;    }    else    {      flux_right = (ch[𝑦] [𝑥+1] – ch[𝑦][𝑥])*𝑛[𝑦][𝑥+1]*𝑔(𝑤[𝑦][𝑥+1])*𝑓(𝑤[𝑦][𝑥+1])/delta𝑋;    }  }  res𝑋 = (flux_left + flux_right)/delta𝑋;  flux_up = 0;  flux_down = 0;  if (𝑦 > 0)  {    if ((ch[𝑦] [𝑥] – ch[𝑦−1] [𝑥]) > 0)    {      flux_up = − (ch[𝑦] [𝑥] – ch[𝑦−1][𝑥])*𝑛 [𝑦−1][𝑥]*𝑔(𝑤[𝑦−1][𝑥])*𝑓(𝑤[𝑦−1][𝑥])/delta𝑌;    }    else    {      flux_up = − (ch[𝑦] [𝑥] – ch[𝑦−1][𝑥])*𝑛 [𝑦] [𝑥]*𝑔 (𝑤[𝑦][𝑥])*𝑓(𝑤[𝑦][𝑥])/delta𝑌;    }  }  if (𝑦 < (size − 1))  {    if ((ch[𝑦+1] [𝑥] – ch[𝑦] [𝑥]) > 0)    {      flux_down = (ch[𝑦+1] [𝑥] – ch[𝑦] [𝑥])*𝑛[𝑦] [𝑥]*𝑔(𝑤[𝑦][𝑥])*𝑓(𝑤[𝑦] [𝑥])/delta𝑌;    }    else    {      flux_down = (ch[𝑦+1] [𝑥] – ch[𝑦] [𝑥])*𝑛[𝑦+1] [𝑥]*𝑔(𝑤[𝑦+1] [𝑥])*𝑓(𝑤[𝑦+1][𝑥])/delta𝑌;    }  }  res𝑌 = (flux_up + flux_down)/delta𝑌;  return (res𝑋 + res𝑌); }

Algorithm 1 

In this code, ch represents the discretization of the proinflammatory cytokine, 𝑛 represents the discretization of neutrophils, 𝑤 is the total density of cells in a position of the space, 𝑥 and 𝑦 are the positions in space, and delta𝑋 and delta𝑌 are the spatial discretizations in 𝑥 and 𝑦 directions, respectively. The test made is to define what is the signal of the characteristic speed, where the speed of the movement of 𝑁(𝑦,𝑥) is given by the term ∇CH. This value is then used to choose between two schemes of finite differences: forward or backward.

We decided to implement our own numerical method to solve the systems of PDEs because (a) we have the possibility to parallelize the code and (b) most of the numerical libraries offer few functions that are suitable to our problem. The sequential code was implemented in C.

9. Numerical Experiments

We performed several simulations in order to verify that the model’s results are in agreement with what is described in the literature. Our objective was to reproduce some characteristics of the microabscess such as an accumulation of dead cells and bacteria in the infection site.

The model’s initial conditions and parameters are given in Tables 1 and 2, respectively.

Table 1 

Initial conditions.

Table 2 

Parameters.

In Table 2, parameters marked with * were adjusted to qualitatively reproduce the results obtained in several studies of the immune response to LPS. In the case of the bacteria (results not shown here), we adjust the equation parameters in order to obtain an exponential decrease in time as shown in [39]. The results of the concentration of proinflammatory cytokines in time (results not shown here) are qualitatively similar to those obtained in some experimental works [40–42]. The parameters marked with ** were based on the values given in the references but were adjusted due to the use of distinct units (e.g., from liter to mm³) or to fit in a 25mm² tissue. In this paper, we obtained parameter values for humans whenever they were available. The variables 𝛽 and total in (5) and (4) were defined based on the concentrations of neutrophils and macrophages per liter given in [43] and were adjusted (a) due to the use of distinct units (e.g., from liter to mm³) and (b) to fit in a 25mm² tissue. The variable total represents the maximum density of cells that fits in a discretized area of the tissue.

In the next sections, we will show the results of the simulation performed with the parameters given in Table 2. In this simulation, we considered a 5mm×5mm two-dimensional domain representing a tissue with 25?mm² of area and a simulation time of 5 days. In our model, the exchange between the vascular system (arterioles and vessels) and tissue was assumed to occur only at the points (1,1), (1,4), and (4,2.5), In this point, immune cells (neutrophils and macrophages) that are in the blood stream can enter into the tissue. The communication between blood vessels and tissue is modeled by permeabilities that vary in time and may depend on the concentration of different cells and molecules (in our model, the endothelium permeability of neutrophils and macrophages depends on the concentration of the proinflammatory cytokine).

10. Bacteria

In the case of bacteria (Figure 1), we observe that initially the bacteria diffuses through the tissue causing tissue damage without its presence to be noticed.

Figure 1 

Temporal evolution and spatial distribution of bacteria.

As soon as resting macrophages residents in the tissue recognize the bacteria they start to produce proinflammatory cytokines that will diffuse through the tissue reaching the blood vessel. Once proinflammatory cytokines interact with the endothelial cells an increase in the endothelium permeability occurs allowing neutrophils and monocytes to migrate to the tissue.

The bacteria starts to die a lot due to the presence of huge numbers of neutrophils. However, the immune response can not completely eliminate bacteria due to the formation of the microabscess pattern. In the microabscess, there are bacteria and a huge concentration of dead cells around it (Figures 4, 5, and 6). In this context, the cleaning process realized by macrophages is very important to allow neutrophils to reach bacteria and eliminate them. Macrophages are responsible for phagocyte dead cells that accumulated in the microabscess. The pattern of microabscess could have lasted longer if we had considered the formation of fibrous tissue around the microabscess.

11. Neutrophil

Neutrophils are initially attracted to the tissue by proinflammatory cytokines produced by activated resident macrophages (Figure 2). Once a neutrophil encounters bacteria, it phagocytizes bacteria and starts to produce proinflammatory cytokines that will attract more neutrophils and macrophages. The cytokine gradient will guide the movement of neutrophils and macrophages in the direction of the highest bacteria concentration.

Figure 2 

Temporal evolution and spatial distribution of neutrophil.

After the microabscess formation, the immune system cells lose contact with a high number of bacteria. As these cells tend to move following the cytokine gradient, we can observe an accumulation of them around the microabscess.

After a significant number of macrophages phagocyte dead cells, the neutrophils can encounter the bacteria and phagocyte them from the border to the center of the microabscess area. As a consequence, a reduction in the microabscess area is observed (Figure 1), which indicates that the immune response is succeeding in controlling the infection.

12. Cytokine

The cytokines in Figure 3 are produced primarily by resident macrophages that are the first to recognize the bacteria presence. The cytokines will increase the endothelium permeability allowing neutrophils to migrate to the tissue. The arrived neutrophils will produce even more cytokines that will guide the movement of neutrophils and macrophage cells in the direction of high concentrations of bacteria.

Figure 3 

Temporal evolution and spatial distribution of cytokine.

Figure 4 

Temporal evolution and spatial distribution of apoptotic neutrophil.

Figure 5 

Temporal evolution and spatial distribution of dead bacteria.

Figure 6 

Temporal evolution and spatial distribution of damaged tissue cells.

During the formation of the microabscess and after it, the production of cytokines is higher in the regions where neutrophils and macrophages have contact with the bacteria that is surrounding the microabscess.

13. Apoptotic Neutrophils

In Figure 4, it can be observed that initially the neutrophils that came from the blood vessel closer to the site of infection died in large number than the neutrophils that came from other sites. Then, after microabscess formation, a lot of neutrophils start to die around the entire microabscess. This phenomenon continues until the microabscess disappears.

14. Dead Bacteria

Figure 5 shows that initially more bacteria died near blood vessels. After the formation of the microabscess, the bacteria starts to die around the entire microabscess since this is the area where the immune response is acting.

15. Tissue Cells

In Figure 6, it can be observed that a lot of tissue damage by bacteria during the period the immune system took to mount an effective immune response. The number of dead tissue cells then reduces, because of the phagocytosis realized by macrophages.

Figure 7 shows the evolution of the healthy tissue area destroyed by the bacteria.

Figure 7 

Temporal evolution and spatial distribution of healthy tissue cells.

16. Microabscess Area

Figure 8 shows the microabscess area using a set of level curves. We defined the microabscess as an area where the concentration of bacteria plus concentration of dead bacteria plus the concentrations of damaged tissue and apoptotic neutrophil is higher. These results show that our model was capable to reproduce the formation of the microabscess in agreement with the observed characteristics of a microabscess [7, 14, 49–54].

Figure 8 

Level curves highlighting the microabscess area at day 5 of the immune response.

17. Conclusions and Future Works

In this work, we presented a computational model for the dynamics of representative types of cells and molecules of the HIS during an immune response to a bacteria. Despite the simplifications and limitations of the model, our results showed that we were able to reproduce an initial microabscess formation. The spatial results show a collection of dead tissue cells, dead bacteria, and apoptotic neutrophil in the microabscess region. This distinct pattern of formation can only be reproduced by spatiotemporal models, such as PDEs models.

As future work, we plan to perform a detailed sensitivity analysis of our microabscess formation model. A previous work [55] has given us some hints about the most sensitivity parameter of the model. We also plan to validate our model against experimental data.

We also plan to modify many aspects of the model to make it more real. For example, we plan to consider a more adequate model to represent the structure of the tissue and its constituents. The tissue can be better characterized as a multiphasic porous medium subjected to stress and deformation variations mainly during the inflammatory process. This porous medium would comprise fluids, extracellular matrix, cells, and molecules. We also plan to model the mechanical behaviour of each of these phases and the mechanical interactions between them.

We have interest in developing models for processes such as vasodilation, coagulation, and others and analyse its effects on the mechanical behaviour of immune system cells and the consequences for the immune response.

With the aim to investigate better the formation of the microabscess, we plan to add another features that contribute to this formation such as the effects of extracellular pH on immune response [56]. Acidic pH predominates at inflammatory loci and other sites of immune activity. Investigations on neutrophils demonstrate mainly inhibition of chemotaxis, respiratory activity, and bactericidal capacity at reduced pH. Besides diminished extracellular pH may play a role in suppressing cytokine production and cytotoxic activities by pulmonary macrophages [56].

Besides we plan to add to the model the process of fibrous tissue formation around the microabscess. We plan to investigate what factors determine if the fibrous tissue will be produced or not. The production of fibrous tissue as well as the coagulation process are ways of the immune system to prevent the bacteria to spread throughout the body doing damage with possible serious consequences, for example, SIRS/MOF. In particular, we are interested in modeling the participation of macrophages, fibroblasts, tissue cells, endothelial cells, and many mediators in the process of fibrous tissue formation. An important step in this process is the production of collagen by fibroblasts induced by the cytokine TGF-𝛽 produced by macrophages [57, 58]. The macrophages has many roles in the processes of wound healing and tissue repair. For example, during the coagulation process, macrophages and endothelial cells are responsible for the production of diverse growth factors and chemotcatic substances that attracts and stimulates the proliferation of tissue cells initiating tissue repair [59].

Acknowledgments

The authors would like to thank FAPEMIG, CNPq, CAPES, and UFJF for supporting this study.