Hybrid approach to model the spatial regulation of T cell responses

Background Moving from the molecular and cellular level to a multi-scale systems understanding of immune responses requires the development of novel approaches to integrate knowledge and data from different biological levels into mechanism-based integrative mathematical models. The aim of our study is to present a methodology for a hybrid modelling of immunological processes in their spatial context. Methods A two-level hybrid mathematical model of immune cell migration and interaction integrating cellular and organ levels of regulation for a 2D spatial consideration of idealized secondary lymphoid organs is developed. It considers the population dynamics of antigen-presenting cells, CD4 + and CD8 + T lymphocytes in naive-, proliferation- and differentiated states. Cell division is assumed to be asymmetric and regulated by the extracellular concentration of interleukin-2 (IL-2) and type I interferon (IFN), together controlling the balance between proliferation and differentiation. The cytokine dynamics is described by reaction-diffusion PDEs whereas the intracellular regulation is modelled with a system of ODEs. Results The mathematical model has been developed, calibrated and numerically implemented to study various scenarios in the regulation of T cell immune responses to infection, in particular the change in the diffusion coefficient of type I IFN as compared to IL-2. We have shown that a hybrid modelling approach provides an efficient tool to describe and analyze the interplay between spatio-temporal processes in the emergence of abnormal immune response dynamics. Discussion Virus persistence in humans is often associated with an exhaustion of T lymphocytes. Many factors can contribute to the development of exhaustion. One of them is associated with a shift from a normal clonal expansion pathway to an altered one characterized by an early terminal differentiation of T cells. We propose that an altered T cell differentiation and proliferation sequence can naturally result from a spatial separation of the signaling events delivered via TCR, IL-2 and type I IFN receptors. Indeed, the spatial overlap of the concentration fields of extracellular IL-2 and IFN in lymph nodes changes dynamically due to different migration patterns of APCs and CD4 + T cells secreting them. Conclusions The proposed hybrid mathematical model of the immune response represents a novel analytical tool to examine challenging issues in the spatio-temporal regulation of cell growth and differentiation, in particular the effect of timing and location of activation signals. Electronic supplementary material The online version of this article (doi:10.1186/s12865-017-0205-0) contains supplementary material, which is available to authorized users.


Numerical method
The T cells zone is considered to be a square domain denoted Ω where the reactiondiffusion equations describing extracellular cytokines concentrations are solved. We apply Dirichlet conditions to the four borders of the square domain (Γ). The numerical implementation of the reaction-diffusion equations was done using the finite difference method. Let us consider the boundary value problem (P ) for the extracellular cytokine field I. We numerically solve the problem using the alternating direction implicit method. We briefly recall the procedure to be followed below.
(P ) : where D is the diffusion coefficient, W and σ are the production and degradation factors respectively, φ(x, y) represents the initial condition condition for I, I 0 is a constant prescribed value at the boundaries. We consider the grid (x i , y j , t n ) = (ih, jh, nδt), where h and δt are the space and time steps respectively. We denote i = 1, 2, ..., N x and j = 1, 2, ..., N y . To begin with, we rewrite the problem (P) in the following form: We split the first equation of the problem (P 1 ) into two sub-steps as follows: We solve the first equation for each fixed j to obtain I n+1/2 . Next, we solve the second to obtain I n . Let us consider the first equation: Rearranging the terms we obtain: Therefore, we can write the first equation of the system 1 in the form: for each fixed j, j = 1, 2, ..., N y − 1, we solve numerically: with the boundary conditions I n i=0 = I 0,1 , I n i=N x = I 0,2 . We solve the equation (2) using Thomas algorithm. For that, we write the left boundary condition I i=0 = I 0,1 as follows: where L 1/2 = 0 and K 1/2 = I 0,1 . Then from the equation (2) for i = 1: . We continue this process for i = 2, 3, ..., N x − 1: where L i+1/2 = −ci bi+aiL i−1/2 , K i+1/2 = fi+aiK i−1/2 bi+aiL i−1/2 . We first obtain the coefficients L i+1/2 , K i+1/2 using the formula (2). Next, we find the solution I n+1/2 for the substep n+1/2 by backward sweep using the equation (3). We apply the same procedure on the second equation of the system (1) to compute the next step solution I n .

Numerical implementation
The source code was written in the Object Oriented Programming (POO) form under C++. The considered time and space units in the model are the minute and the domain length (L) respectively. The library wxWidgets was used to implement a user-friendly interface and visualize the simulations in real-time. The CPU time of numerical simulations was 3 -4 hours on a computer with four cores and 6GB of RAM.

Values of parameters
The number of T cells in the computational domain is ∼ 3×10 3 with the proportions of CD4 + and CD8 + T cells being 2 : 1, and the number of APCs ranging from 30 to 300 cells.
For APC, T-cells and infected cells models, the parameters were fitted with experimental data: 1. k 1 The rate of T-cells production and release into the body: 1.8 /hr; 2. k 2 The death rate of T-cells in the body: 0.12 /hr.