Modelling the Innate Immune System

The Human Immune System (HIS) is a complex network composed of specialized cells, tissues, and organs that is responsible for protecting the organism against diseases caused by distinct pathogenic agents, such as viruses, bacteria and other parasites. The first line of defence against pathogenic agents consists of physical barriers of skin and the mucous membranes. If the pathogenic agents breach this first protection barrier, the innate immune system will be ready for recognize and combat them. The innate immune system is therefore responsible for powerful non-specific defences that prevent or limit infections by most pathogenic microorganisms.


Introduction
The Human Immune System (HIS) is a complex network composed of specialized cells, tissues, and organs that is responsible for protecting the organism against diseases caused by distinct pathogenic agents, such as viruses, bacteria and other parasites.The first line of defence against pathogenic agents consists of physical barriers of skin and the mucous membranes.If the pathogenic agents breach this first protection barrier, the innate immune system will be ready for recognize and combat them.The innate immune system is therefore responsible for powerful non-specific defences that prevent or limit infections by most pathogenic microorganisms.
The understanding of the innate system is therefore essential, not only because it is the first line of defence of the body, but also because of its quick response.However, its complexity and the intense interaction among several components, make this task extremely complex.Some of its aspects, however, may be better understood if a computational model is used.Modelling and simulation help to understand large complex processes, in particular processes with strongly coupled influences and time-dependent interactions as they occur in the HIS.Also, in silico simulations have the advantage that much less investment in technology, resources and time is needed compared to in vivo experiments, allowing researchers to test a large number of hypotheses in a short period of time.
A previous work (Pigozzo et al. (2011)) has developed a mathematical and computational model to simulate the immune response to Lipopolysaccharide (LPS) in a microscopic section of a tissue.The LPS endotoxin is a potent immunostimulant that can induce an acute inflammatory response comparable to that of a bacterial infection.A set of Partial Differential Equations (PDEs) were employed to reproduce the spatial and temporal behaviour of antigens (LPS), neutrophils and cytokines during the first phase of the innate response.
Good modelling practices require the evaluation of the confidence in the new proposed model.An important tool used for this purpose is the sensitivity analysis.The sensitivity analysis consists of the study of the impact caused by the variation of input values of a model on the output generated by it.However, this study can be a time consuming task due to the large number of scenarios that must be evaluated.This prohibitive computational cost leads us to develop a parallel version of the sensitivity analysis code using General-purpose Graphics Processing Units (GPGPUs).GPGPUs were chosen because of their ability to process many streams simultaneously.This chapter describes the GPU-based implementation of the sensitivity analysis and also presents some of the sensitivity analysis results.Our experimental results showed that the parallelization was very effective in improving the sensitivity analysis performance, yielding speedups up to 276.
The remainder of this chapter is organized as follows.Section 2 includes the background necessary for understanding this chapter.Section 3 describes the mathematical model implemented.Section 4 describes the implementation of the GPU version of the sensitivity analysis.Section 5 presents some of the results of the sensitivity analysis and the speedup obtained.Section 7 presents related works.Our conclusions and plans of future works are presented in Section 8.

Biological background
The initial response of the body to an acute biological stress, such as a bacterial infection, is an acute inflammatory response (Janeway et al. (2001)).The strategy of the HIS is to keep some resident macrophages on guard in the tissues to look for any signal of infection.When they find such a signal, the macrophages alert the neutrophils that their help is necessary.The cooperation between macrophages and neutrophils is essential to mount an effective defence, because without the macrophages to recruit the neutrophils to the location of infection, the neutrophils would circulate indefinitely in the blood vessels, impairing the control of huge infections.
The LPS endotoxin is a potent immunostimulant that can induce an acute inflammatory response comparable to that of a bacterial infection.After the lyse of the bacteria by the action of cells of the HIS, the LPS can be released in the host, intensifying the inflammatory response and activating some cells of the innate system, such as neutrophils and macrophages.
The LPS can trigger an inflammatory response through the interaction with receptors on the surface of some cells.For example, the macrophages that reside in the tissue recognize a bacterium through the binding of a protein, TLR4, with LPS.The commitment of this receptor activates the macrophage to phagocyte the bacteria, degrading it internally and secreting proteins known as cytokines and chemokines, as well as other molecules.
The inflammation of an infectious tissue has many benefits in the control of the infection.Besides recruiting cells and molecules of innate immunity from blood vessels to the location of the infected tissue, it increases the lymph flux containing microorganisms and cells that carry antigens to the neighbours' lymphoid tissues, where these cells will present the antigens to the lymphocytes and will initiate the adaptive response.Once the adaptive response is activated, the inflammation also recruits the effectors cells of the adaptive HIS to the location of infection.In CUDA, a parallel function is called kernel.A kernel is a function callable from the CPU and executed on the GPU simultaneously by many threads.Each thread is run by a stream processor.They are grouped into blocks of threads or just blocks.A set of blocks of threads form a grid.When the CPU calls a kernel, it must specify how many threads will be created at runtime.The syntax that specifies the number of threads that will be created to execute a kernel is formally known as the execution configuration, and is flexible to support CUDA's hierarchy of threads, blocks of threads, and grids of blocks.Some steps must be followed to use the GPU: first, the device must be initialized.Then, memory must be allocated in the GPU and data transferred to it.The kernel is then called.After the kernel has finished, results must be copied back to the CPU.

Mathematical model
The model proposed in this chapter is based on a set of Partial Differential Equations (PDEs) originally proposed by Pigozzo et al. (2011).In the original work, a set of PDEs describe the dynamics of the immune response to LPS in a microscopic section of tissue.In particular, the interactions among antigens (LPS molecules), neutrophils and cytokines were modelled.In this chapter, a simplified model of the innate immune system using ODEs is presented to simulate the temporal behaviour of LPS, neutrophils, macrophages and cytokines during the first phase of the immune response.The main differences between our model and the original one (Pigozzo et al. (2011)) are: a) the current model does not consider the spatial dynamics of the cells and molecules and b) the macrophages in two stages of readiness, resting and activated, are introduced in the current model.
Figure 1 presents schematically the relationship between macrophages, neutrophils, proinflammatory cytokines and LPS.LPS cause a response in both macrophages and neutrophils, that recognize LPS and phagocyte them.The process of phagocytosis induces, in a rapid way, the apoptosis of neutrophils.This induction is associated with the generation of reactive oxygen species (ROS) (Zhang et al. (2003)).The resting macrophages become activated when they find LPS in the tissue.The pro-inflammatory cytokine is produced by both active macrophages and neutrophils after they recognize LPS.It induces an increase in the endothelial permeability allowing more neutrophils to leave the blood vessels and enter the infected tissue.
Our set of equations is given below, where RM, AM, A, N and CH represent the population of resting macrophages, activated macrophages, LPS, neutrophils and pro-inflammatory cytokines, respectively.The dynamics of LPS is modelled with Equation 1.
The term µ A A models the decay of LPS, where µ A is its decay rate.The term −(λ N|A .N + λ AM|A .AM + λ RM|A .RM).A models the phagocytosis of LPS by macrophages and neutrophils, where λ N|A is the phagocytosis rate of neutrophils, λ AM|A is the phagocytosis rate of active macrophages, and λ RM|A is the phagocytosis rate of resting macrophages.
Neutrophils are modelled with Equation 2.
The term permeability N uses a Hill equation (Goutelle et al. (2008)) to model how permeability of the endothelium of the blood vessels depends on the local concentration of cytokines.Hill equations are also used, for example, to model drug dose-response relationships (Wagner (1968)).
The idea is to model the increase in the permeability of the endothelium according to the concentration of the pro-inflammatory cytokines into the endothelium.In the Hill equation, P max N represents the maximum rate of increase of endothelium permeability to neutrophils induced by pro-inflammatory cytokines, P min N represents the minimum rate of increase of endothelium permeability induced by pro-inflammatory cytokines and keqch is the concentration of the pro-inflammatory cytokine that exerts 50% of the maximum effect in the increase of the permeability.The term µ N N models the neutrophil apoptosis, where µ N is the rate of apoptosis.The term λ A|N A.N models the neutrophil apoptosis induced by the phagocytosis, where λ A|N represent the rate of this induced apoptosis.The term source N represents the source term of neutrophil, that is, the number of neutrophils that is entering the tissue from the blood vessels.This number depends on the endothelium permeability (permeability N ) and the capacity of the tissue to support the entrance of neutrophils (N max ), that can also represent the blood concentration of Neutrophils.
The dynamics of cytokine is presented in Equation 3.
The term µ CH CH models the pro-inflammatory cytokine decay, where µ CH is the decay rate.
The term (β CH|N N + β CH|AM AM).A models the production of the pro-inflammatory cytokine by the neutrophils and activated macrophages, where β CH|N and β CH|AM are the rate of this production by neutrophils and macrophages, respectively.
Equation 4 presents the dynamics of the resting macrophages.
The term permeability RM models how permeability of the endothelium of the blood vessels to macrophages depends on the local concentration of cytokines.The term µ RM RM models the resting macrophage apoptosis, where µ RM is the rate of apoptosis.
Finally, the dynamics of activate macrophages is presented in Equation 5.
The term µ AM RM models the activated macrophage apoptosis, where µ RM is the rate of apoptosis.

Implementation
The sensitivity analysis consists in the analysis of impacts caused by variations of parameters and initial conditions of the mathematical model against its dependent variables (Saltelli et al. (2008)).If a parameter causes a drastic change in the output of the problem, after suffering a minor change in its initial value, it is thought that this parameter is sensitive to the problem studied.Otherwise, this variable has little impact in the model.The sensitivity analysis is used to improve the understanding of the mathematical model as it allows us to identify input parameters that are more relevant for the model, i.e. the values of these parameters should be carefully estimated.In this chapter we use a brute force approach to exam the influence of the 19 parameters present in the equation and two of the initial conditions.A small change in the value of each parameter is done, and then the model is solved again for this new parameter set.This process is done many times, since all combinations of distinct values of parameters and initial conditions must be considered.We analyse the impact of changing one coefficient at a time.The parameters and initial conditions were adjusted from -100% to + 100% (in steps

355
Modelling the Innate Immune System www.intechopen.com of 2%) of their initial values, except for some parameters, that were also adjusted from -100% to + 100%, but in steps of 20%.The combination of all different set of parameters and initial conditions give us a total of 450,000 system of ODEs that must be evaluated in this work.
The sequential code that implements the sensitivity analysis was first implemented in C.
Then the code was parallelized using CUDA.The parallel code is based on the idea that each combination of distinct values of parameters and initial conditions can be computed independently by a distinct CUDA thread.The number of threads that will be used during computation depends on the GPU characteristics.In particular, the number of blocks and threads per block are chosen taking into account two distinct values defined by the hardware: a) the warp size and b) the maximum number of threads per block.
The forward Euler method was used for the numerical solution of the systems of ODEs with a time-step of 0.0001 days.The models were simulated to represent a total period equivalent to 5 days after the initial infection.

Experimental evaluation
In this section the experimental results obtained by the execution of both versions of our simulator of the innate system, sequential and parallel, are presented.The results reveal that our CUDA version was responsible for a significant improvement in performance: a speedup of 276 was obtained.This expressive gain was due to the embarrassingly parallel nature of computation that must be performed.In particular, the same computation must be performed for a huge amount of data, and there are no dependency and/or communication between parallel tasks.

Simulation
To study the importance of some cells, molecules and processes in the dynamics of the innate immune response, a set of simulations were performed for distinct values of parameters and initial conditions.The complete set of equations that has been simulated, including the initial values used, are presented by Equation 6: It should be noticed that in this case two distinct initial values for A(0) will be used: A(0)= 20 and A(0)=40.
The sensitivity analysis has shown that two parameters are relevant to the model: the capacity of the tissue to support the entrance of new neutrophils (N max ) and the phagocytosis rate of LPS by neutrophils (λ N|A ).
N max is the most sensitive parameter in the model.The capacity of the tissue to support the entrance of new neutrophils is directed related to the permeability of the endothelial cells, which form the linings of the blood vessels.If a positive adjustment is made in the parameter related to the permeability, then there are more neutrophils entering into the tissue.This larger amount of neutrophils into the tissue has many consequences: first, more cells are phagocyting, so the amount of LPS reduces faster.Second, a smaller amount of resting macrophages becomes active, because there is less LPS into the tissue.Third, a larger amount of cytokines are produced, since neutrophils are the main responsible for this production.If a negative adjustment is made, the inverse effect can be observed: with a smaller amount of neutrophils in the tissue, more resting macrophages become active.Also, a smaller amount of cytokines are produced.In the second scenario, with the double of LPS and starting with just one resting macrophage, it can be observed that bringing more neutrophils into the tissue do not reduce the number of resting macrophages that become active.This happens due to the larger amount of LPS in this scenario when compared to the previous one.The larger amount of activated macrophages also explains why the amount of cytokines in this scenario is larger than in the previous one.Figures 7 to 11 present the complete scenario.
The third scenario presents the results obtained when the initial amount of LPS is again equal to 20.This scenario revels that the second most sensitive parameter is λ N|A .λ N|A is responsible for determining how effective is the phagocitosis of the neutrophils in tissue.It can be observed in Figures 12 to 16 that a negative adjustment in this tax makes the neutrophil response to be less effective against LPS, while a positive adjustment in the tax makes the neutrophil response to be more effective.Resting macrophages and activated macrophages are also affected by distinct values of λ N|A .Increasing the value of λ N|A causes the neutrophils  to produced more cytokines, so more macrophages can migrate into the tissue through blood vessel, and also there are more cells into the tissue that can phagocyte LPS.
The last scenario is presented by Figures 17 to 21.In this scenario, the amount of LPS is doubled when compared to the previous one.It can be observed that distinct values used as initial conditions for LPS only changes how long it takes to the complete elimination of LPS.It can also be observed that both macrophages populations are affected by the larger amount of LPS.In particular, the amount of macrophages is slightly higher in this scenario due to the larger amount of LPS.

Related works
This section presents some models and simulators of the HIS found in the literature.Basically two distinct approaches are used: ODEs and PDEs.

ODEs models
A model of inflammation composed by ODEs in a three-dimensional domain considering three types of cells/molecules has been proposed by Kumar et al. (2004): the pathogen and two inflammatory mediators.The model was able to reproduce some experimental results depending on the values used for initial conditions and parameters.The authors described the results of the sensitivity analysis and some therapeutic strategies were suggested from this analysis.The work was then extended (Reynolds et al. (2006)) to investigate the advantages of an anti-inflammatory response dependent on time.In this extension, the mathematical model was built from simpler models, called reduced models.The mathematical model (Reynolds et al. (2006)) consists of a system of ODEs with four equations to model: a) the pathogen; b) the active phagocytes; c) tissue damage; and d) anti-inflammatory mediators.
A new adaptation of the first model (Kumar et al. (2004)) was proposed to simulate many scenarios involving repeated doses of endotoxin (Day et al. (2006)).In this work the results

361
Modelling the Innate Immune System www.intechopen.comA one-dimensional model to show if and when leukocytes successfully defend the body against a bacterial infection is presented in Keener & Sneyd (1998).A phase-plane method is then used to study the influence of two parameters, the enhanced leukocyte emigration from bloodstream and the chemotactic response of the leukocytes to the attractant.
Finally, one last work (Vodovotz et al. (2006)) developed a more complete system of ODEs of acute inflammation, including macrophages, neutrophils, dendritic cells, Th1 cells, the blood pressure, tissue trauma, effector elements such as iNOS, NO − 2 and NO − 3 , pro-inflammatory and anti-inflammatory cytokines, and coagulation factors.The model has proven to be useful in simulating the inflammatory response induced in mice by endotoxin, trauma and surgery or surgical bleeding, being able to predict to some extent the levels of TNF, IL-10, IL-6 and reactive products of NO (NO − 2 and NO − 3 ).2011) present a PDE model to simulate the immune response to lipopolysaccharide (LPS) in a microscopic section of a tissue, reproducing, for this purpose, the initiation, maintenance and resolution of immune response.

Other works
Several proposals which attempt to model both the innate and the adaptive HIS can be found in the literature.An ODE model is used to describe the interaction of HIV and tuberculosis with the immune system (Denise & Kirschner (1999)).Other work focus on models of HIV and T-lymphocyte dynamics, and includes more limited discussions of hepatitis C virus (HCV), hepatitis B virus (HBV), cytomegalovirus (CMV) and lymphocytic choriomeningitis virus (LCMV) dynamics and interactions with the immune system (Perelson (2002)).An ODE model of cell-free viral spread of HIV in a compartment was proposed by Perelson et al. (1993).
Another interesting work tries to integrate the immune system in the general physiology of the host and considers the interaction between the immune and neuroendocrine system  ( Muraille et al. (1996)).Klein (1980) presents and compares three mathematical models of B cell differentiation and proliferation.
ImmSim (Bezzi et al. (1997); Celada & Seiden (1992)) is a simulator of the HIS that implements the following mechanisms: immunological memory, affinity maturation, effects of hypermutation, autoimmune response, among others.CAFISS (a Complex Adaptive Framework for Immune System Simulation) (Tay & Jhavar (2005)) is a framework used for modelling the immune system, particularly HIV attack.SIMMUNE (Meier-Schellersheim & Mack (1999)) allows users to model cell biological systems based on data that describes cellular behaviour on distinct scales.Although it was developed to simulate immunological phenomena, it can be used in distinct domains.A similar tool is CyCells (Warrender (2004)), designed to study intercellular relationships.

Conclusion and future works
In this chapter we presented the sensitivity analysis of a mathematical model that simulates the immune response to LPS in a microscopic section of a tissue.The results have shown that the two most relevant parameters of the model are: the capacity of the tissue to support the entrance of more neutrophils and the phagocytosis rate of LPS by neutrophils.
The sensitivity analysis can be a time consuming task due to the large number of scenarios that must be evaluated.This prohibitive computational cost leads us to develop a parallel version of the sensitivity analysis code using GPGPUs.Our experimental results showed that the parallelization was very effective in improving the sensitivity analysis performance, yielding speedups up to 276.

Fig. 1 .
Fig. 1.Relationship between the components.(GPGPUs).CUDA includes C software development tools and libraries to hide the GPGPU hardware from programmers.
Initial conditions, parameters and units.

Fig. 2 .
Fig. 2. Temporal evolution of cytokines with A(0)=20 and for distinct values of N max .

359
Modelling the Innate Immune System www.intechopen.com

Fig. 5 .
Fig. 5. Temporal evolution of resting macrophages with A(0)=20 and for distinct values of N max .

Fig. 7 .
Fig. 7. Temporal evolution of cytokines with A(0)=40 and for distinct values of N max .

Fig. 9 .
Fig. 9. Temporal evolution of LPS with A(0)=40 and for distinct values of N max .
Computational Algorithms and Their Applications www.intechopen.com

Fig. 11 .
Fig. 11.Temporal evolution of activate macrophages with A(0)=40 and for distinct values of N max .
Computational Algorithms and Their Applications www.intechopen.com

Fig. 15 .
Fig. 15.Temporal evolution of resting macrophages with A(0)=20 and for distinct values of λ N|A .

Table 1 .
The experiments were performed on a 2.8 GHz Intel Core i7-860 processor, with 8 GB RAM, 32 KB L1 data cache, 8 MB L2 cache with a NVIDIA GeForce 285 GTX.The system runs a 64-bits version of Linux kernel 2.6.31 and version 3.0 of CUDA toolkit.The gcc version 4.4.2 was used to compile all versions of our code.The NVIDIA GeForce 285 GTX has 240 stream processors, 30 multiprocessors, each one with 16KB of shared memory, and 1GB of global memory.The number of threads per block are equal to 879, and each block has 512 threads.The codes were executed 3 times to all versions of our simulator, and the average execution time for each version of the code is presented in Table1.The standard deviation obtained was negligible.The execution times were used to calculate the speedup factor.The speedup were obtained by dividing the sequential execution time of the simulator by its parallel version.Serial and parallel execution times.All times are in seconds.
Table 2 presents the initial conditions and the values of the parameters used in the simulations of all cases.Exceptions to the values presented in Table 2 are highlighted in the text.
Figures 2 to 6 illustrate this situation.It can be observed that the LPS decays faster when N max achieves its maximum value.