Edited by: Gennady Bocharov, Institute of Numerical Mathematics (RAS), Russia
Reviewed by: Christian Kurts, University of Bonn, Germany; Raffaele De Palma, Università degli Studi della Campania Luigi Vanvitelli Caserta, Italy
This article was submitted to T Cell Biology, a section of the journal Frontiers in Immunology
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
The adaptive immune response is initiated in lymph nodes by contact between antigen-bearing dendritic cells (DCs) and antigen-specific T cells. A selected number of naïve T cells that recognize a specific antigen may proliferate into expanded clones, differentiate, and acquire an effector phenotype. Despite growing experimental knowledge, certain mechanistic aspects of T cell behavior in lymph nodes remain poorly understood. Computational modeling approaches may help in addressing such gaps. Here we introduce an agent-based model describing T cell movements and their interactions with DCs, leading to activation and expansion of cognate T cell clones, in a two-dimensional representation of the lymph node paracortex. The primary objective was to test the putative role of T cell chemotaxis toward DCs, and quantitatively assess the impact of chemotaxis with respect to T cell priming efficacy. Firstly, we evaluated whether chemotaxis of naïve T cells toward a nearest DC may accelerate the scanning process, by quantifying, through simulations, the number of unique T cell—DC contact events. We demonstrate that, in the presence of naïve T cell-to-DC chemoattraction, a higher total number of contacts occurs, as compared to a T cell random walk scenario. However, the forming swarm of naïve T cells, as these cells get attracted to the neighborhood of a DC, may then physically restrict access of additional T cells to the DC, leading to an actual decrease in the cumulative number of unique contacts between naïve T cells and DCs. Secondly, we investigated the potential role of chemotaxis in maintaining cognate T cell clone expansion. The time course of cognate T cells number in the system was used as a quantitative characteristic of the expansion. Model-based simulations indicate that inclusion of chemotaxis, which is selective for already activated (but not naïve) antigen-specific T cells, may strongly accelerate the time of immune response occurrence, which subsequently increases the overall amplitude of the T cell clone expansion process.
After maturation in the thymus, immunologically-naïve T lymphocytes (or T cells) continuously circulate between the blood and secondary lymphoid organs, including lymphatic nodes (LNs) and the spleen. Each one of the millions of T cell clones bear unique T cell receptors (TCRs), which define their antigen specificity. In LNs, naïve T cells may encounter dendritic cells (DCs) presenting cognate antigens as MHC-bound peptides (pMHC) on their surface. As a result of such a specific and durable contact T cell-to-DC contact, a naïve T cell may become activated and subsequently proliferate and differentiate into effector forms. This constitutes, in most simplified terms, the essence of the immune response. The fraction of naïve T cells that recognize a particular antigen can be as small as 10−5-10−6 (
Over the past two decades, experiments using two-photon microscopy (2PM) have been applied in the study of murine LNs
Intravital LN observations have shown that cognate T cell interactions with antigen-presenting DCs can be categorized into several stages—with may possibly overlap over time (
Despite such detailed observations, there is no comprehensive understanding, yet, of the detailed mechanisms and dynamics of immune cell interactions; in particular, the fate of individual cells is difficult to track for longer periods of time
Overall, an agent-based model (ABM) represents a system of interest, with a definition of key players and of relevant interactions among these players that influence the system's behavior. A typical ABM consists of a simulation space (world), stand-alone objects (agents), and rules to set the behavior of individual agents as well as interactions among them (rules). Thus, in our 2D ABM framework, we consider T cell motility and emerging interactions of immune cells within the LN T zone. The T zone was modeled as a lattice of 100 × 100 patches, resulting in an effective physical surface area of 500 × 500 μm2 (i.e., 5 × 5 μm2 per patch). This modeling framework also considers two types of agents: T cells and antigen-presenting DCs.
In order to reproduce interactions between agents present in the LN T zone, 2,000 naïve T cells were randomly placed in this square domain, along with 8 antigen-bearing DCs, randomly placed in 8 fixed positions, as shown in
2D ABM model simulation snapshots. T cells are featured as small filled circles of different colors. DCs are shown as larger black cross symbols. “Chemokine cloud” areas are colored in gray, around DCs. Green bars on the top and bottom boundaries of the computational domain physically restrict T cell movement (no escape), while yellow spaces, in between, represent MS structures (T cell escape allowed).
For simplification, antigen-presenting DCs were assumed to be immobilized, given their migration speed in LNs vs. that of T cells is relatively slow (
The stochastic motion of a T cell was implemented in two ways, depending on the particular mechanism tested in the model:
To capture the short-term persistence character of T cell movement, T cell motion at each time step was set to be dependent on its previous direction. The new direction was thus calculated by defining a take-off angle from the previous direction, as the sum of two random angles sampled clockwise and counterclockwise from a uniform distribution in the range of (0, θmax) degrees. If the resulting adjacent lattice position happened to be already occupied (e.g., by another T cell, or a DC, or a boundary patch on the top or bottom side of the computational domain), then the T cell was computed to remain in its current position, and another direction calculation would be attempted at the next time step.
In our simulations, we aimed at reproducing not only realistic velocities of T cell movement, but also a physiologically-relevant T cell motility coefficient—realistically reflecting the short-term persistence character of T cell movement. Motility coefficients used in such 2D simulations were calculated from time lapse microscopy records using the formula:
where 〈Δ
After some time spent in the “chemokine cloud,” a T cell loses its sensitivity to the chemokine gradient and starts moving randomly again. The time of T cell de-sensitization to chemokine(s) was taken to be 10 min (
Periodic boundary conditions were applied on the left and right sides of the computational domain: if a T cell were to leave the domain through one side, it would be allowed to immediately re-appear from the other side of the domain, moving in the same direction. In contrast, T cells were not allowed to randomly cross top and bottom boundaries of the computational domain. These boundaries, instead, contained “open patches,” functionally corresponding to medullary sinuses (MS) and efferent lymphatics in a LN. Accordingly, if a T cell were to leave the computational domain through an MS patch at either the top or bottom side, a new T cell would be allowed to enter the computational domain, from a random position through its lateral borders. These settings allowed us to keep an overall constant T cell density in the system under study. The model was explored using a range of numbers of patches containing such MS escape structures (8–120 patches).
Relatively rare cognate naïve T cells, capable of recognizing DC-presented antigens, were included in the model. No difference in movements, between cognate vs. non-cognate naïve T cells before their first encounter with DC was assumed. After a first encounter, a cognate T cell was set to form a stable contact with a DC for ~24 h; specifically, the duration of contact for each particular cognate T cell was a random value generated from a normal distribution with a mean of 24 h and a variance of 2 h. Upon completion of such a contact, a T cell became activated. The model subsequently simulated activated T cells which randomly migrated into a virtual LN and interacted with DCs bearing cognate antigen complexes. We also explored the option of chemoattraction of activated (but not naïve) T cells toward a neighboring DC.
Similarly to the models by Bogle et al. (
where
For activated T cells which ended up outside a DC contact zone, the stimulation level was set to decrease according to an exponential law:
where λ is the exponent indicator, corresponding to a half-life period of 24 h (
Thus, a T cell was set to divide when the stimulation level
Two factors were used in the model, to limit the proliferation of activated T cells: (i) a minimal time of about 8 h (random value generated from a normal distribution with a variance of 1 h for each newly formed T cell) was set between successive T cell divisions; (ii) a maximal number of successive proliferating T cell divisions was set. Following the last division of an activated T cell reaching effector status, no further division was allowed, and the cell was eliminated from the system within 24 h. Also, to avoid a strong increase in overall T cell density in the model (as a result of cognate T cell expansion), the following rule was added: if an activated T cell were to leave the computational domain
The following quantitative outcome measures were generated:
LN transit time: average time from the moment a T cell object appears in the computational domain of the model, until it exits that space through MS patches.
Total number of T cell–DC contacts: a cumulative count of contact events between any DC-T cell pair, including possible repeat contacts, as detected during a given simulation time.
Calculation of the number of unique T cell–DC contacts: only the first contact of a given T cell with any of the DCs was taken into account. In the model, all DCs were assumed to be strictly identical.
Dynamics of cognate T cell number in the computational domain: included naïve, activated and effector T cells, in simulations of up to Day 28.
Cumulative outflux of cognate T cells: this was taken as the characteristic measure of the adaptive immune response intensity.
To calculate prediction intervals (90% PI) for each of the model outcomes, 45–100 independent ABM simulations were performed using identical model parameter values, yet different, randomly generated initial T cell and DC positions within the computational domain.
The
Using the 2D ABM model of a LN as described above, a large number of simulations were performed for various model parameter settings. For simplification purposes in these exploratory simulations, several model parameters were fixed using biologically reasonable estimates. In particular, the assignment of values to T cell and DC densities (respectively, 2,000 and 8 cells per 500 × 500 μm2) and to the T cell movement parameter reflecting short-term persistence (θmax = 80°) allowed us to reproduce,
Our first goal was to explore the impact of chemoattraction upon efficiency in the process of T cell repertoire scanning, as naïve T cells moved toward a DC. The measure of such efficacy was computed as the rate of accumulation of unique T cell—DC contact events. The chemotaxis strength estimate (
We specifically explored the model behavior in relationship to the parameter reflective of the size of overall medullary sinuses (MS) and corresponding to the number of “open patches” at the top and bottom boundaries of the computational domain. This parameter is critical in the model, as it regulates T cell turnover rate and thus influences most of the quantitative outcomes. As shown in
As shown in
Dependence of the number of total
As shown in
A second goal of this study was to explore simulations of cognate T cell priming and expansion under different model parameter settings. Here, we used small, albeit non-zero values for cognate T cell frequency parameters, i.e., the probability of a new incoming T cell to be cognate, which could be activated in contact with DCs and proliferate as described above. Also, we tested different maximal numbers of activated T cell divisions, i.e., maximal numbers of T cell generations starting from naïve cognate T cells to finally differentiated effector T cell. Proceeding from the previous part, we fixed the overall MS size at 32 patches: it ensured a transit time through the LN T zone in agreement with experimental values. All simulations were performed under the assumption of a constant level of antigen stimulation in the LN: neither the number of DCs, nor their positions, nor their properties changed during the simulations.
In preliminary simulations not shown here, we sought to reproduce T cell immune responses using the same “random walk” scenario for both naïve and activated lymphocytes; under such conditions, cognate T cell expansion levels were low and not robust. By visual inspection of such simulation trajectories, we observed that most activated cognate T cells would leave the computational domain prior to any cell division occurring. To resolve this technical modeling issue, we allowed chemotaxis toward a neighboring DC to be selective for already activated (but not naïve) cognate T cells. As shown in
Simulations of cognate T cell numbers in the LN T-zone
In particular, following days 3–5 of simulations (
Model outcomes were highly sensitive to the activated T cell chemotaxis strength value (
Another important factor was the cognate T cell frequency, i.e., the probability for new incoming naïve T cells to recognize antigens presented by DCs (
Only moderate increases in “steady-state” cognate T cell numbers and in the corresponding outflows were observed, when the parameter value reflecting the maximal number of T cell divisions was increased from 10 to 20 (
The original motivation driving this modeling study was to enable the exploration of dynamic spatial effects, in particular a detailed investigation of the relationship between T cell motility behavior and the timing and intensity of an immune response. Intravital microscopy (2PM) yields a wealth of information within a very restricted region of the LN and during a short period of time (hours), whereas histology provides complementary views, yet limited to two dimensions and with no dynamic time element. An ABM of the LN allows for integrative simulations which may help filling the gaps between these two experimental approaches. To explore several hypotheses on T cell motility and their interactions with DCs, we independently developed a simplified, yet biologically reasonable version of a 2D ABM of the LN T cell zone. The model was developed and qualified using
The 2D ABM presented here, which includes T cell movement, activation and proliferation in a LN, allows for the integration of a number of processes and at the scale of the LN system, and illustrates the necessity to consider all essential processes simultaneously, in order to generate a realistic dynamic picture of the immune response. Because detailed experimental data required to characterize all these processes are not available, we made assumptions regarding several parameters embedded in the model; nonetheless, the model developed here realistically reproduced key temporal characteristics, such as the T cell motility coefficient (
Operating characteristics of the model were supported by sensitivity analyses, whereby simulations were run over ranges of model parameter values, with model outcomes being compared against physiological values. Model simulations indicated that the generated T cell response was sensitive to factors, such as naïve cognate T cell frequency and the strength of the hypothetical chemoattraction of T cells toward neighboring DCs. Thus, despite a simplified, semi-empirical structure of the model, we obtained reasonable and robust simulations over a wide range of unknown parameter values.
The ABM LN model presented here was set as a two-dimensional (2D) model, rather than a more physiological three-dimensional (3D) model. A 2D setting of the model allowed us, obviously, to drastically reduce the computational cost of ABM simulations, while carefully estimating the impact of stochastic effects on simulation outcomes. In most of ABMs, T cell movement is implemented as a “random walk” in a non-guided biophysical domain; under such settings, advantages of 3D vs. 2D models may, in fact, not be obvious. The role of the fibroblastic reticular cell (FRC) network in guiding T cell motion in the LN has been studied over many years (
The numerical simulations reported here were focused on two pivotal questions. The first question focused on whether local chemoattraction of T cells toward DCs would promote or hamper the scanning efficiency of DCs, within a LN. We demonstrated that, for an effective DC scanning of the T cell repertoire, a T cell “random walk” motility scenario appeared to be the optimal strategy (
Such results are in full agreement with results from earlier 2D ABM research, in which T cell motility was accurately captured to help determine the impact of chemotactic attraction of T cells-to-DC on repertoire scanning (
The relevance of a chemo-attraction process on T cell scanning efficiency by DCs was also addressed in modeling work by Vroomans et al. (
The second question which we sought to address here was about a potential role for chemotaxis in immune response initiation. For such a purpose, we simplified the description of cognate T cell activation, to minimize the overall number of parameters in the model. The concept of a TCR stimulation signal (
Using our 2D ABM approach, we determined that a feature of selection for activated cognate T cells is required, to reproduce their pronounced expansion upon response to antigen stimulation. As mentioned above, activated T cell chemoattraction lead to T cell accumulation in swarms, around DCs. This effectively caused a longer half-life of these T cells in the LN, slowing down their elimination from MS patches, and also causing frequent contacts with DCs, thereby contributing to activation signaling (
Similarly to previous ABM applications tailored to the LN, we showed that such a modeling technique proves to be a useful tool to integrate current knowledge and data on molecular and cellular interactions between immune cells, to then generate novel hypotheses which may guide further experimental studies, to overall improve our mechanistic understanding of the immune activation process that takes place in the LN. Many questions on this dynamical process in the LN remain open, in particular questions on the emerging role of the FRC network in regulating immune responses. Future developments of 3D models, with detailed stromal elements, may play an important role in further elucidating biological mechanisms (
IA and YK designed, initiated and developed the modeling project. GH and KP participated in discussions and reviewed modeling design and concept. YK and IA wrote the computer program and translated the biological model into agent-based modeling language. IA performed model calculations. GH elaborated on paper outline and points of focus. All authors contributed toward manuscript writing and revisions.
The authors declare that this study received funding from AstraZeneca. IA, KP, YK are employed by and KP, YK are owners of M&S Decisions, a modeling consultancy which received research funding from AstraZeneca. GH is employed by, and shareholder of AstraZeneca.
We are most grateful to Dr. Uri Wilensky (Northwestern University) for distributing the
The Supplementary Material for this article can be found online at:
2D agent-based modelling movie at the start of simulation.
2D agent-based modelling movie at the 7th day post numerical experiment start.
Additional details on the model development and analysis, e.g. table with parameter values, results of sensitivity analysis and additional model simulations.
NetLogo 5.0.1 scripts corresponding to the final model version presented in the manuscript.