# Spatial Lymphocyte Dynamics in Lymph Nodes Predicts the Cytotoxic T Cell Frequency Needed for HIV Infection Control

Dmitry Grebennikov, Anass Bouchnita, Vitaly Volpert, Nikolay Bessonov, Andreas Meyerhans, Gennady Bocharov

<figure><img src="image_1.png" /></figure>

## To cite this version:

Dmitry Grebennikov, Anass Bouchnita, Vitaly Volpert, Nikolay Bessonov, Andreas Meyerhans, et al. Spatial Lymphocyte Dynamics in Lymph Nodes Predicts the Cytotoxic T Cell Frequency Needed for HIV Infection Control. Frontiers in Immunology, 2019, 10, pp.1213. 10.3389/fimmu.2019.01213. [hal-02363399](https://hal.science/hal-02363399v1)

HAL Id: hal-02363399

https://hal.science/hal-02363399v1

Submitted on 13 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Spatial Lymphocyte Dynamics in Lymph Nodes Predicts the Cytotoxic T Cell Frequency Needed for HIV Infection Control

Dmitry Grebennikov<sup>1,2,3*</sup>, Anass Bouchnita<sup>4</sup>, Vitaly Volpert<sup>3,5,6</sup>, Nikolay Bessonov<sup>7</sup>, Andreas Meyerhans<sup>8,9</sup> and Gennady Bocharov<sup>2,10*</sup>

**OPEN ACCESS**

**Edited by:**
Michael Loran Dustin,
University of Oxford, United Kingdom

**Reviewed by:**
Jens Volker Stein,
Université de Fribourg, Switzerland
Martin Meier-Schellersheim,
National Institutes of Health (NIH),
United States
Frederik Graw,
Universität Heidelberg, Germany

*Correspondence:
Dmitry Grebennikov
dmitry.ew@gmail.com
Gennady Bocharov
bocharov@m.inm.ras.ru

**Specialty section:**
This article was submitted to
*T Cell Biology,*
a section of the journal
*Frontiers in Immunology*

**Received:** 08 February 2019

**Accepted:** 13 May 2019

**Published:** 11 June 2019

**Citation:**
Grebennikov D, Bouchnita A, Volpert V, Bessonov N, Meyerhans A and Bocharov G (2019) Spatial Lymphocyte Dynamics in Lymph Nodes Predicts the Cytotoxic T Cell Frequency Needed for HIV Infection Control. Front. Immunol. 10:1213. doi: 10.3389/fimmu.2019.01213

The surveillance of host body tissues by immune cells is central for mediating their defense function. *In vivo* imaging technologies have been used to quantitatively characterize target cell scanning and migration of lymphocytes within lymph nodes (LNs). The translation of these quantitative insights into a predictive understanding of immune system functioning in response to various perturbations critically depends on computational tools linking the individual immune cell properties with the emergent behavior of the immune system. By choosing the Newtonian second law for the governing equations, we developed a broadly applicable mathematical model linking individual and coordinated T-cell behaviors. The spatial cell dynamics is described by a superposition of autonomous locomotion, intercellular interaction, and viscous damping processes. The model is calibrated using *in vivo* data on T-cell motility metrics in LNs such as the translational speeds, turning angle speeds, and meandering indices. The model is applied to predict the impact of T-cell motility on protection against HIV infection, i.e., to estimate the threshold frequency of HIV-specific cytotoxic T cells (CTLs) that is required to detect productively infected cells before the release of viral particles starts. With this, it provides guidance for HIV vaccine studies allowing for the migration of cells in fibrotic LNs.

**Keywords:** lymphoid tissue, cell motility, HIV infection, cytotoxic T cell scanning, multicellular dynamics, dissipative particle dynamics, stochastic differential equation

## INTRODUCTION

The surveillance of host body tissues by cells of the immune system is central for mediating defense functions against invading pathogens and tumor cells (1, 2). The initial recognition of foreign antigens that leads to the induction of adaptive immune responses takes place in lymph nodes (LNs), which, by virtue of their location and structure, facilitate the interactions between immune cells (3). The motility of pathogen spread and immune cells represents relevant parameters

<footnote n="1"><sup>1</sup> Moscow Institute of Physics and Technology, National Research University, Dolgoprudny, Russia, <sup>2</sup> Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences, Moscow, Russia, <sup>3</sup> Peoples' Friendship University of Russia (RUDN University), Moscow, Russia, <sup>4</sup> Division of Scientific Computing, Department of Information Technology, Uppsala University, Uppsala, Sweden, <sup>5</sup> Institut Camille Jordan, UMR 5208 CNRS, University Lyon 1, Villeurbanne, France, <sup>6</sup> INRIA Team Dracula, INRIA Lyon La Doua, Villeurbanne, France, <sup>7</sup> Institute of Problems of Mechanical Engineering, Russian Academy of Sciences, Saint Petersburg, Russia, <sup>8</sup> Infection Biology Laboratory, Department of Experimental and Health Sciences, Universitat Pompeu Fabra, Barcelona, Spain, <sup>9</sup> Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Spain, <sup>10</sup> Sechenov First Moscow State Medical University, Moscow, Russia</footnote>
controlling the fate of the pathogen-host interaction. *In vivo* imaging technologies have been used to quantitatively characterize target cell scanning and migration dynamics of lymphocytes within LNs (4, 5). The translation of these quantitative insights into a predictive understanding of immune system functioning in response to various perturbations critically depends on the availability of computational tools linking the individual immune cell properties with the systems response as a whole (6).

Multiscale models of the immune system provide the *in silico* tool to embed immune processes into their spatial context (7-9). A core module of the models is the mathematical framework used to describe individual cell migration in complex multicellular environments. One can distinguish two general types of modeling approaches, cellular automata-based models (CAMs), and physical models (PMs). CAMs consider a regular grid with cells that change their state in time and space according to some rules (functions of the system state). The respective computational algorithms can take the form of random walks (10) or cellular Potts models (11). Although CAMs incorporate experimentally defined characteristics of cell motion and, thus, simulate cell dynamics based on actual data, they lack quantifiable links to the underlying biophysical interactions between cells in multicellular environments and to intrinsic cell motility parameters (12). PMs of lymphocyte migration dynamics derived from the Newtonian second law offer the possibility to define cell motions in terms of the forces generated by the environment and the cell itself. Using the experimental data on cell movement, the potential functions underlying cell-to-cell interactions and intrinsic cell motility can be identified and can provide a deeper insight into the mechanical properties of cells. Thus, PMs of individual cells and coordinated cell migration represent a general and generic way to describe and predict the multicellular system dynamics for a broad range of cell numbers and external conditions (13, 14).

It is widely accepted in immunology that the physiological function of cytotoxic T cell (CTL) motility is to search for target cells, i.e., for virus-infected cells or cancer cells (15). Computational modeling studies have revealed that the search efficiency depends on the organization of the stromal environment of a tissue (16). In addition, the spatial behavior, for example, of HIV-infected target cells scanned for foreign antigens by CTLs strongly impacts the elimination efficiency of the infected targets (17, 18). Experimental investigation of live attenuated SIV vaccines clearly suggested that a robust protection against intravenous wild-type SIVmac239 challenge strongly correlates with the number and function of antigen-specific effector CTLs in LN rather than the responses of such cells in the blood (19). However, the quantitative effects of T-cell migration parameters in LNs on the efficiency of antiviral immune responses *in vivo* remain unknown.

In the current study, we have developed a physics-based description of spatial T-lymphocyte dynamics in the multicellular environment of LNs. A fundamental relationship between a cell motion and the forces acting on it is provided by Newton's second law. It is used to formulate, calibrate, and apply a generic mathematical model of coordinated T-cell migration

dynamics in LNs. By choosing a first principles approach in formulating the governing equations in conjunction with published experimental data on T-cell motility in lymphoid tissues, we offer a broadly applicable generic mathematical tool linking individual and coordinated cell behaviors. The potential of the model is illustrated by an analysis of the combined effects of antigen-specific T-cell numbers and intrinsic T-cell motility parameters in LNs on the time needed to locate both mobile and non-motile HIV-infected target cells. Computed predictions of the ratio of effector CTLs to infected T cells in the LN paracortex needed for a timely detection of infected cells within 18 h postinfection, i.e., before the release of viral particles starts (20), provide a novel quantitative guide for an informed design of HIV vaccines.

# MATERIALS AND METHODS

## Programming Languages and Computing Resources

All algorithms were written in C++ and compiled using G++ (version 5.4.0). Pseudorandom numbers were generated using the PCG random library (version 0.98) and the PCG64-XSL-RR algorithm (21). The seed was either specified manually (for code development) or set based on the system's random device (for computational experiments). Simulations were run on a 2-core Xeon E3-1220 v5 @3.0 GHz × 4 processor. The wxWidgets library (version 2.8.12) was used for visualization purposes. The processing of the simulation results (i.e., calculating statistical motility profiles, comparing CDFs, and plotting) was implemented in Python and R scripts.

## Model Equations of Multicellular Dynamics

According to a basic mechanics view, a system consisting of N cells of some mass located in a liquid milieu, interacting with each other and affected by some external field, is uniquely determined by their coordinates and velocities and is governed by the classical mechanics motion equations. In our model, each cell $i, i = 1, \bar{N}$, is represented as the circle with certain mass $m_i$, radius $r_i$, and position of its center $x_i$. The fundamental equation governing locomotion of cells is Newton's second law of motion. It can be expressed as follows:

$$m_i \ddot{x}_i = F_i = \sum_{j \neq i} f_{ij}^{\text{int}} + f_i^{\text{mot}} + f_i^{\text{dis}}, \quad i = \overline{1, \bar{N}}, \quad (1)$$

where the first term on the right side specifies the net effect of the pairwise interaction forces with contacting neighbor cells, the second term stands for the cell intrinsic locomotion force, by which the cell establishes motility within the extracellular matrix (ECM) of the LN reticular network, and the last one takes into account the action of a dissipative force, taken to be proportional to the cell velocity $f_i^{\text{dis}} = -\mu \dot{x}_i$. We neglect the impact of gravity.

## Random Motility Force Sampling

The random motility force $f_i^{\text{mot}}$ for the ith cell is modeled as a stochastic vector $f_i$ sampled every 30 s from certain probability distributions analogously to the inverse homogeneous correlated
random walk (IHomoCRW) model (22). The motility magnitude $|f_i| = \eta_i \cdot K_i$ is sampled from the following Gaussian distribution: $K_i \in |N(\mu(K), \sigma^2(K))|$. To obtain the motility magnitude $|f_i|$, the sampled value $K_i$ is multiplied by the arresting coefficient $\eta_i$. The arresting coefficients are increased for both T cells and DCs if they establish a sufficiently long contact to temporarily arrest their inner motility as follows: (1) $\eta_i := 10\eta_i^{\text{default}}$ for T cells and DCs when the duration of an uninterrupted contact exceeded 30 s, and (2) $\eta_i := 100\eta_i^{\text{default}}$ if the contact outlasted 20 min. The cell inner motility is restored back to a default value if the contact lasted for a time longer than the sampled value $t_{\text{contact}} \in N(2, 0.4)$ hours. The parameter $\eta_i$ is also used to decrease intrinsic motility when performing in silico simulations to study the effect of decreased T-cell motility on target cell location efficiency (see details in **Supplementary Text**).

The motility direction $\hat{f}_i$ is turned from the previous direction on the angle $\theta_i$:

$$ \alpha_i \in N(0, \sigma^2(\alpha)), \quad \theta_i = \alpha_i \cdot \left(1 - \left(\frac{K_i}{K_{\text{max}}}\right)^{\beta}\right), \quad K_{\text{max}} = \mu(K) + 3\sigma(K). \quad (2) $$

Here, $N(0, \sigma^2(\cdot))$ denotes a Gaussian distribution, and $\beta$ is a scalar coefficient. The angle sampled from the normal distribution is multiplied by a factor depending on the sampled motility magnitude to reproduce the experimentally observed negative correlation between cell translational and turning angle speeds. Indeed, the cells do not simultaneously perform fast translational movements and large reorientations (22). Note that a similar feature was named “directional propensity” and modeled with trigonometric parameterization in a cellular Potts model to describe the motion of T cells (11). The Gaussian distribution for the motility magnitude is set so that the $(\mu - 3\sigma, \mu + 3\sigma)$ range is positive. The absolute value is taken to ensure that the magnitude is non-negative. The parameter $K_{\text{max}}$ provides an upper boundary for sampled values $K_i$ (approximately 1 of 370 cases falls outside of the three-sigma interval). The hat above the vector denotes the normalized unit vector.

## Implementation of Contact Inhibition of Locomotion

After the stochastic vector $f_i$ is sampled, it is modified in accordance with the contact inhibition of locomotion (CIL) model, as described (23). The resultant vector $f_i^{\text{mot}}$ is then used in the right-hand side of Equation (1). The modification consists of shifting the direction of vector $f_i$ away from the neighboring cells and decreasing the magnitude of vector $f_i$ proportionally to the number of neighboring cells:

$$ f_i^{\text{mot}} = \frac{|f_i| \cdot (c_{\text{inh}}\hat{f}_i + \hat{R}_i)}{c_{\text{inh}} + n}, \quad \hat{R}_i = \sum_{j, h_{ij} \le r_i + r_j} \frac{x_i - x_j}{h_{ij}}, \quad (3) $$

in which $|f_i|$ is the magnitude and $\hat{f}_i$ is the direction of the inner motility as it would be if unaffected by CIL, $n$ is the number of

neighboring cells in contact (such that the distance between cell centers $h_{ij} \le r_i + r_j$), and $\hat{R}_i$ determines the net shift of the inner motility direction away from the neighboring cells, $c_{\text{inh}}$ is the weighting coefficient varying the level of CIL. The hat above the vector indicates that it is normalized.

## Numerical Integration of the Equations of Cell Motion

To numerically integrate the equations of motion (Equation 1), we used the first-order semi-implicit (i.e., the cell coordinate at time $t^{n+1}$ is computed using the velocity vector $v_i^{n+1}$ rather than $v_i^n$) Euler method:

$$ v_i^{n+1} = \frac{m_i v_i^n + h \cdot (F_i^{\text{int}}(t^n, x_i^n) + f_i^{\text{mot}}(t^n, x_i^n))}{m_i + h \cdot \mu} \quad (4) $$

$$ x_i^{n+1} = x_i^n + h \cdot v_i^{n+1} \quad (5) $$

in which $x_i^n$ and $v_i^n$ are the coordinate and velocity of cell $i$ at the time $t^n$ after $n$ steps $t^n = t^0 + h \cdot n$. We note that the second-order generalization of this method, i.e., the Störmer–Verlet method, could be developed. However, it will be computationally more demanding as the cell acceleration depending on velocity due to the presence of dissipative velocity-damping viscosity forces needs to be reevaluated at each time step $t^{n+1}$. We verified that the time step $h = 0.02$ min used in the simulations is sufficient for a stable integration of the initial value problem with the semi-implicit Euler method. To efficiently locate the neighboring cells (which is needed for intercellular force calculations and for determining the effect of CIL), we use a simple uniform-grid-based spatial neighbor search, which performs well for a densely packed multicellular environment. Note that the convergence of the integration scheme was verified by repeating simulations for a smaller time step.

## Boundary Conditions

During the model calibration process, we used periodic boundary conditions for all boundaries of a square domain. To perform in silico simulations in a closed ellipse-shaped domain representing a LN, we implemented a biologically based boundary condition of cell repolarization. We do not model explicitly the interaction forces between cells and the boundary (i.e., the subcapsular sinus wall). At the stage of coordinate updates (in accordance with the numerical scheme specified in the section Numerical Integration of the Equations of Cell Motion), if the proposed coordinate of cell $x_i^{n+1}$ is outside the boundary, the current coordinate of the cell is preserved ($x_i^{n+1} = x_i^n$), while the direction of the motility vector $\hat{f}_i^{\text{mot}}$ is changed to be the opposite direction of vector $v_i^{n+1}$, thus resulting in cell repolarization.

## Generating the Initial Spatial Configuration for Simulations Within a LN

To generate the initial spatial configuration of the immune cells within a LN, we followed the descriptions from a LN imaging study (24). The following cell subsets were considered: CD4$^+$ T cells, CD8$^+$ T cells, and cross-presenting migratory CD8$\alpha^{\text{int}}$CD103$^{\text{hi}}$ DCs. Both T-cell subsets are distributed
uniformly through the whole LN, while migratory DCs are found mainly deep in the paracortex area. To arrange cells in agreement with the experimental data, we approximated the DC-rich area as an ellipse $\Omega_{DC}^{\alpha=0.99}$. The spatial positions for DC locations are iteratively sampled from the 2D Gaussian distribution with a 99-percentile ellipse $\Omega_{DC}^{\alpha=0.99}$ and accepted if the DC with sampled coordinates lies within the LN domain $\Omega_{LN}$ and does not overlap with the other seeded DCs. After DCs are placed, the T cells are positioned uniformly through the remaining non-occupied space of $\Omega_{LN}$.

# RESULTS

## Biophysical Parametrization of the Spatial Multicellular Dynamics

Multicellular systems dynamics can be accurately described by biophysical models as reviewed recently (14, 25). Here, we develop a physics-based mathematical model of coordinated immune cell motion that belongs to the class of self-propelled particle models (14) and, more generally, to the dissipative particle dynamics (26, 27) framework.

Immune cells in LNs are continuously interacting with each other and with stromal cells via forces of different origin, i.e., elastic (membranes), chemical (receptors), and electric. The respective forces in combination with cell intrinsic locomotion events act in concert to determine the basal intranodal motility of T cells. Figure 1A presents the overall summary of physical forces included in the model with some implementation details. The scheme of the forces exerted on cell i interacting with cells p and k is shown in Figure 1B. The quantitative features of the force functions are detailed in Figure 1C. Here, $f_{ij}^{\text{int}}$ is the intercellular force acting on cell i due to interaction with cell j. The pairwise cell-to-cell interactions are assumed to have a finite cutoff distance and are considered to be elastic acting along the line of cell centers. The intercellular forces $f_{ij}^{\text{int}}$ can be considered as the gradients of pairwise potentials, which are repulsive at short distances and attractive at larger distances, thus accounting for volume exclusion at the cell body and cell-cell adhesion near membranes. We consider the following cubic polynomial function to model the force exerted by cell j on cell i:

$$f_{ij}^{\text{int}} = \frac{x_i - x_j}{h_{ij}} \cdot \begin{cases} -a \cdot f^{\text{adh}} \cdot \frac{r_j - x}{r_j} + b \cdot f^{\text{adh}} \cdot \left( \frac{r_j - x}{r_j} \right)^3, & h_{ij} < r_i + r_j, \\ 0, & h_{ij} \ge r_i + r_j, \end{cases} \quad (6)$$

where $r_i$ is a radius of the ith cell membrane, $h_{ij}$ is the distance between cell centers (see Figure 1C), and $x = h_{ij} - r_i$ is the distance between the center of cell j and the surface membrane of cell i. The function $a \cdot f^{\text{adh}} \cdot \frac{r_j - x}{r_j}$ describes the attraction force between two cells, and the function $b \cdot f^{\text{adh}} \cdot \left( \frac{r_j - x}{r_j} \right)^3$ corresponds to a repulsive force, both calibrated as shown in Figure 1C. The coefficients a and b are set such that the minimum of function $f_{ij}^{\text{int}}$ is equal to $f^{\text{adh}}$. Thus, the only remaining free parameter is $f^{\text{adh}}$, the adhesive interaction strength. In the case of T cell/T cell interaction, it corresponds to weak non-specific electrical forces (electrostatic and electrodynamic) that

are expected to be present between all cells according to the model of Bell (28). We calibrate this parameter by the typical value of low-adhesive forces, with which integrins present on T-cell membrane bind to their ligands present on the other cells (29). For cognate T cell/APC interactions the attraction force is much stronger as it is determined by a broad spectrum of various adhesion molecules involved in T-cell activation clusters, i.e., the immunological synapse (30). The estimated values of the intercellular interaction forces are given in Table 1. For details on the data-based T-cell motility model calibration, see Supplementary Text.

The dissipative (friction) force acting on T cells describes the effect of viscous damping, which reduces the velocity of the cell. It is assumed to be proportional to the cell velocity $f_i^{\text{dis}} = -\mu \dot{x}_i$. The dissipative force acts along the line of the cell center and in opposite direction to the cell displacement. Consideration of viscous damping is appropriate for the highly viscous low-Reynolds-number environment of LNs (40). The viscous damping parameter estimate is listed in Table 1.

The random motility force $f_i^{\text{mot}}$ determines the traction of self-propelled lymphocytes. It represents a stochastic process of receptor-mediated cell-ECM interactions regulated by either cytoskeletal or membrane reorganizations and governed by biomechanical and intracellular molecular mechanisms (4, 13). Basically, cells establish directed caterpillar-like movement by polarizing, forming contacts between their leading edge and collagen fibers of ECM, detaching their trailing edge from ECM, and contracting. However, T lymphocytes and dendritic cells (DCs) are characterized by low-adhesive integrin interactions with the microenvironment. This allows them to adapt their direction and morphology with no need to reorganize microstructure while effectively sliding along the stromal network of fibroblastic reticular cells (41, 42). As we do not model the reticular network and the ECM microstructure in this study explicitly, this motility behavior is considered implicitly in the stochastic nature of $f_i^{\text{mot}}$. Note that the autonomous cell motility can also be affected by external signaling, e.g., through chemotaxis, CIL, or immunological synapse formation. The cell trajectory in the model is characterized by three quantifiable values, i.e., the translational speed, the turning angle speed, and the meandering index as explained in Figure 1D and as described in Read et al. (22). The corresponding experimental data are shown in Figure 1E. To capture the experimentally observed patterns of T-lymphocyte migration in lymphoid tissues (see Figure 1E), the T-cell motility is modeled using a random variable $f_i$ with its magnitude and angle values updated every $\Delta t$ seconds according to the IHomoCRW recently suggested and validated (22). The IHomoCRW model was shown to reproduce the experimentally measured statistical profiles of T-cell locomotion (22). In the present model, the magnitude and direction of the random vector $f_i^{\text{mot}}$ are sampled from distributions provided by the experimental data (the specific rules are defined in section Materials and Methods). The key difference from the original IHomoCRW model is that it is the cell motility inducing force $f_i^{\text{mot}}$ rather than the cell velocity $\dot{x}_i$ that is sampled and then substituted into equation (Equation 1). In addition,
A

<table><thead><tr><td>Force</td><td>Description</td><td>Implementation</td></tr></thead><tbody><tr><td>Self-propelling</td><td>Stochastic process of active intrinsic lymphocyte motility</td><td>Sampling the motility force vector from probability distributions according to lHomoCRW model</td></tr><tr><td>Intercellular</td><td>Attractive cell-cell adhesion and repulsive volume exclusion + contact inhibition of locomotion</td><td>Calibrated pairwise forces between contacting cells + decreasing motility force and shifting it from neighbor cell aggregates</td></tr><tr><td>Dissipative</td><td>Friction with viscous LN environment</td><td>Viscous damping proportional to the cell velocity</td></tr></tbody></table>

B

$$m_i \ddot{x}_i = F_i^{\text{int}} + f_i^{\text{mot}} + f_i^{\text{dis}} = \sum_{j \neq i} f_{ij}^{\text{int}} + f_i^{\text{mot}} - \mu \ddot{x}_i$$

D

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

C

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

G

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

E

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

Protocol: * 30-min pre-run to randomly mix cells
* 10 experiments of 30 minutes length
* save cell positions every 30 seconds
* exclude cells with total displacement < 27 µm
* exclude cells passed through periodic boundaries

<caption>**FIGURE 1** | Physics-based model of multicellular system dynamics reproduces experimental data on T-cell locomotion. (A) The set of forces considered in the model with description of their features and implementation details. (B) The fundamental equation governing locomotion of cells determined by the forces exerted on cell *i*, (Continued)</caption>
**FIGURE 1** | including the repulsive–attractive interaction with neighbor cells *p* and *k*, respectively. (C) The parameterization of intercellular interaction force $f_{ij}^{int}$ and formula definition. The calibrated force for non-specific interaction of two T cells with a radius of 3 µm is depicted. By simulation, the parameters *a* and *b* are calculated at each time step depending on the radii $r_i, r_j$ and the distances $h_{ij}, x$, so that the condition $f_{ij}(\lambda r_j) = f_{ij}(r_j) = 0, \min_{ij}(x) = -f_{ij}^{adh}$ is satisfied. The parameter $\lambda$ determines the relative deformation of the cells that separates the repulsive and attractive interactions between them. Parameter $f_{ij}^{adh}$ represents the adhesive strength between the membranes of cells *i* and *j*. (D) Schematic illustration and definition of the metrics characterizing T-cell motility: translational speed, turning angle speed, and meandering index. All metrics are measured for each cell every $\Delta t$ seconds and pooled together to form statistical distributions. (E) Statistical profiles characterizing the T-cell locomotion consists of distribution histograms of translational speeds, turning angle speeds, and meandering indices. The histograms are derived from the corresponding empirical cumulative distribution functions (CDFs) available in Figure S17 from Read et al. (22), in which original *in vivo* data are presented. (F) The details of the 2D geometric setup for simulations used in the model calibration: spatial configuration, initial and boundary conditions, and the experimental protocol used to sample the statistical profile. (G) The statistical characteristics of T-cell motility coming from simulations of the calibrated model plotted against the *in vivo* histogram data (22). The statistical distributions of each metric are depicted as CDFs. The Kolmogorov–Smirnov statistics comparing the model and target CDFs are indicated with their respective *p*-values.

<table><caption>**TABLE 1** | Set of calibrated model parameters used as a baseline for all simulations.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>References</th></tr></thead><tbody><tr><td>m<sub>TC</sub></td><td>T-cell mass</td><td>215 pg</td><td>(31–34)†*</td></tr><tr><td>r<sub>TC</sub></td><td>T-cell radius</td><td>3 µm</td><td>(35)†</td></tr><tr><td>μ</td><td>Viscous damping coefficient</td><td>0.2 nN·min/µm (= 12 g/s)</td><td>(36, 37)†</td></tr><tr><td>f<sub>T–T</sub><sup>adh</sup></td><td>Adhesive strength between T cells</td><td>0.01 nN</td><td>(28–30)†</td></tr><tr><td>λ<sub>TC</sub></td><td>The normalized distance between the cell-to-cell interaction synapse and the cell center separating repulsive and attractive modes of T-cell interaction (model analog of nuclear-to-cytoplasmic ratio in experimental cell biology)</td><td>0.83 (with respect to r<sub>TC</sub>)</td><td>(38)†</td></tr><tr><td>Δt</td><td>Time step for inner motility m<sub>TC</sub> update</td><td>30 s</td><td>(22)‡</td></tr><tr><td>μ (m<sub>TC</sub>)</td><td>Mean of the inner motility force magnitude distribution</td><td>3 nN</td><td>(22, 39)‡†</td></tr><tr><td>σ (m<sub>TC</sub>)</td><td>SD of the inner motility force magnitude distribution</td><td>0.3 nN</td><td>(22, 39)‡†</td></tr><tr><td>σ (α<sub>TC</sub>)</td><td>SD of the inner motility turning angle distribution</td><td>60°</td><td>(22)‡‡</td></tr><tr><td>c<sub>inh</sub></td><td>Scalar coefficient varying the level of CIL</td><td>1.0</td><td>(23)‡‡</td></tr><tr><td>β</td><td>Scalar coefficient in parameterization of negative correlation between magnitude and turning angle of sampled inner motility</td><td>2.0</td><td>(23)‡‡</td></tr></tbody></table>

the random vector $f_i^{mot}$ can be influenced by contact effects from neighboring cells, resulting in (1) a shift of the vector $f_i^{mot}$ away from neighboring cells and (2) a decrease of its magnitude proportionally to their number, similar to the CIL model (23) (see details in section Materials and Methods). By default, the arresting coefficient for T cells is equal to one. For DCs, its value is estimated so that the resultant DC velocities do not exceed 5 µm/min (the estimated value is specified in Table 2).

Overall, the mechanistic model of the spatial multicellular dynamics is formulated as a system of N random ordinary differential equations (44) represented by Equation (1) and embedded into the 2D geometric domain as detailed in Figure 1F. Essentially, the system is a deterministic system of ordinary differential equations on each interval of $\Delta t$ seconds, until the force $f_i^{mot}$ becomes updated. The quantitative consistency of the computational model of multicellular dynamics with experimental data on translation speed, turning angle speed, and the meandering index is illustrated in Figure 1G. The relevant components of the numerical implementation of the model (computational domain, boundary conditions, integration algorithm) are described in Materials and Methods. The dynamics of the net forces and their contributions

acting on a randomly selected T cell in a simulation of multicellular dynamics are shown in Figure S1.

## Calibration of T-Cell Motility

Our model mostly operates with biophysical parameters that are either directly measurable or can be estimated indirectly such as the mass *m* (wet weight) and the radius *r* of a cell, the adhesive strength between T-cell membranes $f_{ij}^{adh}$ (measured by single cell force spectroscopy), the viscous damping coefficient $\mu$, typical forces and velocities of T cells, and the location of demarcation between repulsive and attractive areas of a cell $\lambda$ (nuclear-to-cytoplasm ratio). The other parameters that describe the random motility force or the contact inhibition of locomotion are derived using the information presented in the original IHomoCRW model (22) and the CIL model (23) with the underlying experimental data. To calibrate our model, we evaluated admissible ranges of parameters and tuned them manually to match the statistical characteristics of T-cell locomotion (22). The baseline sets of the estimated model parameters are presented in Tables 1, 2. For details of the parameter estimation, see Supplementary Text.

<footnote>†Parameters obtained directly from experimental measurements.</footnote>

<footnote>*Parameters estimated indirectly from experimental measurements.</footnote>

<footnote>‡Parameters derived from underlying computational models.</footnote>

<footnote>‡Parameters tuned to fit cell motility profiles within the model calibration.</footnote>
<table><caption>**TABLE 2** | Extra parameters which are used for LN simulations.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value/range</th><th>References</th></tr></thead><tbody><tr><td>m<sub>CD4</sub></td><td>CD4<sup>+</sup> T-cell mass</td><td>215 ± 28 pg</td><td>(31–34)</td></tr><tr><td>m<sub>CD8</sub></td><td>CD8<sup>+</sup> T-cell mass</td><td>290 ± 28 pg</td><td>(31–34)</td></tr><tr><td>m<sub>DC</sub></td><td>Dendritic cell mass</td><td>350 ± 28 pg</td><td>(32–34)</td></tr><tr><td>r<sub>DC</sub></td><td>Dendritic cell radius</td><td>6.5 µm</td><td>(38, 43)</td></tr><tr><td>λ<sub>DC</sub></td><td>The normalized distance between the cell-to-cell interaction synapse and the cell center separating repulsive and attractive modes of DC interaction (model analog of nuclear-to-cytoplasmic ratio in experimental cell biology)</td><td>0.5</td><td>(38)</td></tr><tr><td>f<sub>adh</sub><sub>DC</sub></td><td>Adhesive strength between specific T cells and DCs</td><td>1 nN</td><td>(28, 30, 36)</td></tr><tr><td>η<sub>DC</sub></td><td>Default value for coefficient arresting inner motility of DCs</td><td>3</td><td>Tuned so that DC velocities are &lt;5 µm/min</td></tr></tbody></table>

## Computational Domain, Immune Cell Subsets, and Initial Configuration

The computational domain was implemented as an ellipse-shaped 2D approximation of the bean-like cross section of a murine skin-draining LN (see **Figure 2A**). At the beginning, both CD4⁺ T cells (green) and CD8⁺ T cells (blue) are evenly distributed throughout the domain. Some randomly chosen T cells are considered to be antigen specific and marked in light green and blue, respectively (their numbers are specified below). The antigen-presenting cells considered in this study represent the subset of cross-presenting migratory CD8αintCD103hi DCs, which are mainly involved in CD8⁺ T-cell immune responses and which immigrate into LNs from the periphery (24). They are normally localized in the deep parts of the T-cell zone and leave LNs slowly with a turnover rate of 6 days. For initial configuration, these DCs are spatially placed according to a Gaussian distribution with 99-percentile ellipse $\Omega_{DC}^{\alpha=0.99}$ representing the T-cell zone (see **Figure 2A** and section Materials and Methods).

The numbers of antigen-specific DC and T-cell subsets are estimated using published data (24), which were rescaled according to the size of the computational domain. A total population of 12,469 immune cells was considered. The total number of non-antigen-specific T cells was estimated so that about 80% of the computation domain was filled up. The precursor frequency of antigen-specific T cells, that is, their proportion in the total amount of T cells, was set to be about 1%. We consider the inflow and outflow of immune cells to the region of interest to be negligible because of the short simulation time of 12 h. The closed boundary conditions used in the simulations are specified in section Materials and Methods. The overall geometrical scheme of the computational domain and the initial configuration of the multicellular system generated for simulations are presented in **Figure 2A**.

## Data Assimilation and Model Validation

To assimilate the statistical data on the three T-cell locomotion measures (i.e., the translational speed, turning angle speed, and meandering index), the following numerical simulation protocol was used, which is close to the original experimental protocol (22). First, the same 2D 412 × 412 µm² domain was used, in which we initialized 4,489 squarely tiled T cells with 3-µm radii and $\eta \approx 80\%$ packing density. The initial direction of the intrinsic

motility force was generated randomly for all cells. The positions of cells were saved every 30 s during 10 numerical experiments of 30-min simulation time after a 30-min pre-run to randomly mix the cells. Cells with total displacements <27 µm were excluded as was done in the original experimental protocol. Likewise, cells that passed through the boundary and left the imaging volume were also excluded. The saved cell positions were post-processed to calculate the target metrics (defined in **Figure 1D**), which were pooled together to form three separate distributions. The pooled cell motility distributions were calibrated with *in vivo* data. The simultaneous adjustment of all distributions was computationally challenging due to the different uncorrelated aspects of cell migration captured in each of the motility metric as previously outlined (22). **Figure 1G** shows the best-fit cumulative distribution functions (CDFs) of the calibrated model with the baseline parameter set from **Table 1** and the target experimentally observed distributions with Kolmogorov-Smirnov statistics and *p*-values describing the discrepancy between CDFs.

The evolution of the above multicellular system was simulated over a 12-h period. The visualization of the systems spatiotemporal dynamics is presented in **Movie S1**. **Figures 2B,C** shows the kinetics of median velocities of antigen-specific CD4⁺ T and CD8⁺ T cells, and the median distances between the T cells and the centroid of their cognate antigen-presenting DCs throughout 12 h of an *in silico* experiment. The model demonstrates that antigen-specific CD8⁺ T cells that interact with their cognate CD8αint DCs but not the CD4⁺ T cells decrease their velocities, move closer to the area of DCs in the first 4–6 h, and remain there with low velocities afterward. **Figures 2B,C** is quantitatively consistent with experimental data shown in **Figures 1E,F**, and in Figure 2B from Kitano et al. (24).

## Quantitation of the DC and T-Cell Contact Interactions

The calibrated mathematical model of T-cell locomotion was validated by confronting its predictions with data from the intranodal spatiotemporal dynamics of different immune cell subsets after soluble antigen immunization presented in a recent experimental study (24). The data specify the evolution of the distances between the centroid of the migratory DC area and individual CD4⁺ T and CD8⁺ T cells. The model was adjusted to the functional configuration of skin-draining LNs specified
<figure><img src="image_14.png" /></figure>

in the above study. A representative example of the numerical simulation of individual cell trajectories is shown in **Figure 3A**. An example of multicellular dynamics in a LN during 12 h is shown in **Figure 3B**.

## Quantitation of the Forces Determining T-Cell and DC Motility and Their Interaction

To consider DCs in multicellular system simulations, we carried out a parameterization of their intrinsic motility forces and the intercellular forces for contacts between (1) two DCs, (2)

antigen-presenting DCs and antigen-specific T cells, and (3) antigen-presenting DCs and polyclonal T cells. The values of the corresponding parameters are presented in **Table 2**. The physical forces driving the dynamics of individual cells in the LN and the respective velocities of the cells predicted by the model are shown in **Figure 3C**. We assume that the intrinsic motility of DCs can be represented by the same type of force function as that for T cells (**Figure 1C**); however, due to their much smaller average velocity, the respective DC force function value was of small magnitude. The adhesive force for cognate contacts
<figure><img src="image_15.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>D</figcaption></figure>

<figure><img src="image_20.png" /></figure>

Quantity and duration of cognate DC - CD8+ TC contacts

<figure><img src="image_21.png" /></figure>

**FIGURE 3** | Quantitation of immune cell motility, driving forces, and contacts. (A) Representative example of individual cell trajectories obtained with numerical realization of the calibrated model. The trajectories illustrate the 5-h dynamics of 15 cells randomly chosen from 4,489 cells presented in the 412 × 412 µm² domain with periodic boundary conditions. (B) Twelve-hour multicellular dynamics of T-cell trajectories in a lymph node obtained by numerical simulation with an initial configuration specified in Figure 2A. Only cells with total displacement longer than 27 µm are shown. (C) Values of forces and cell velocities driving the multicellular system dynamics in a square subdomain of a LN. In a center pane, the velocity field is represented as a contour plot of the field of cell velocity magnitudes linearly interpolated at uniform grid, as well as detected streamlines of possible cell flow patterns. (D) Kinetics of the numbers of cognate DC-T cell contacts at different stages of the simulation and distribution of durations of all cognate contact durations occurring within a 12-h simulation. DC, dendritic cell.

(i.e., of antigen-specific T cells with antigen-presenting DC) is around 100 times higher (~1 nN) than the non-specific adhesion force for T cell/T cell contacts (36). We also implemented a

computational procedure to temporarily arrest the motility for T cells in a sufficiently long cognate contact (see section Materials and Methods).
## Ag-Specific CD8+ T Cells Migrate Toward Cross-Presenting DCs and Form Cognate Contacts With Them

**Figure 3D** presents the model prediction for the kinetics of the number of cognate DC–CD8+ T-cell contacts occurring at different time intervals during the *in silico* simulation. Antigen-specific CD8+ T cells robustly increase the number of contacts with DCs over time in the process of T-cell zone scanning for antigen-presenting target cells. Although most of the cognate contacts are of short duration, i.e., they last for <5 min, the distribution has a heavy tail of stable more than 1-h length contacts. These predictions are in agreement with previous data (45).

## CTL Frequency Needed to Locate HIV-Infected Target Cells Before Viral Release

During viral infections, the induction of cellular immune responses takes place in secondary lymphoid organs such as LNs and spleen. Antigen-presenting cells such as DCs take up antigens and migrate to LNs to encounter specific lymphocytes, e.g., CD4+ T and CD8+ T cells, to induce their activation and differentiation into effector and memory cell subtypes (46). The low frequency of antigen-specific T cells in unprimed hosts turns the scanning of cognate DCs by specific T cells in a highly crowded cellular LN environment into a “needle-in-a-haystack” problem (47). It was revealed that optimal LN scanning depends on a combination of intrinsic T-cell motility, the chemokine milieu, and the microarchitecture of the LN (1). When virus-infected DCs reach the LN, the less the time needed to locate virus-specific T cells and to form stable DC–T cell contacts, the more likely is that the precursor CTL activation will happen before the viruses will be released from infected cells, therefore making the elimination of local clusters of infection spread more probable. This aspect of CD8+ T-cell activity is crucial for a prompt activation of specific CTL immune responses and the elimination of viruses. The precursor frequency in blood can be as small as 0.0001% (48), reaching about 5–10% in the chronic stage of an HIV infection (49). The here-developed physics-based model of T-cell dynamics can be directly used to study the efficiency of scanning the paracortical T-cell zone of the LN for target cells expressing cognate antigen as a function of the frequency of CTL and their motility.

Development of an effective AIDS vaccine remains a global priority, and there is a need for a vaccine to induce cellular immune responses capable of eradicating or efficiently containing virus replication (50). Experimental studies with attenuated SIV vaccines indicated that SIV-specific CTLs, if present in sufficient frequencies, can completely control and even clear an infection (19). Similar to SIV, HIV infection is sustained by the activation of CD4+ T cells, which occurs in the form of transient bursts in the local microenvironment of lymphoid tissues (51, 52). The proximal activation and transmission involving latently infected cells represent locally propagating events (53). Therefore, we applied our calibrated model of spatial immune cell dynamics in LNs to study the necessary conditions for effector HIV-specific CTLs to promptly locate HIV-infected

target cells before they can release viral progeny. We consider only one HIV-infected cell in the computational domain, which is consistent with the frequency of productively infected CD4+ T cells of about 0.0001–0.001 (54). Specifically, the newly infected target cell should be located by the nearby effector cells before it can release viral progeny, i.e., before completion of the 18–24 h life cycle of HIV (20).

The overall simulation setup is the same as described in the model validation subsection above. Randomly chosen cells in the stochastically generated multicellular system configurations representing the LN cortex zone were marked as infected in yellow (see **Figure 4A**). Both the motile CD4+ T cells and the non-motile DCs were considered as HIV-infected targets. In simulations, we varied the frequency of HIV-specific CD8+ T cells and the intrinsic motility of T cells (searching CD8+ T cells, infected- and uninfected CD4+ T cells) (**Figure 4A**) to analyze the effect of variations on the target cell detection time. A 10-fold range of HIV-specific CD8+ T-cell frequencies typical for HIV infection, i.e., from 0.4 to 5%, was examined. The intrinsic motility of T cell was varied within 100 and 50% relative to the calibrated baseline parameters of average T-cell velocity (see details in the **Supplementary Text**). A decreased intranodal T-cell motility (of searching CD8+ T cells, infected- and uninfected CD4+ T cells) is expected to take place during the chronic stage of an HIV infection when LN tissues become fibrotic, i.e., when collagen formation in T-cell zones takes place (55). Then, T cells have to move through increased collagen deposition with major consequences for search patterns (56). In addition, CD4+ T-cell migration is also inhibited by the HIV-1 Nef protein as shown in chemotaxis assays (57). In our study, the motility of all considered types of T-cell subsets, i.e., the searching CD8+ T cells, uninfected CD4+ T cells, and infected CD4+ T cells, is decreased uniformly.

**Figures 4B,C** illustrate the model predictions for the decrease of time to locate HIV-infected target cells with the increase of HIV-specific CD8+ T-cell frequency. The modeling results imply that 5% is a sufficient effector CTL frequency for a timely detection of both types of target cells within 18 h post-infection, i.e., before the beginning of HIV release from productively infected cells (20). A stepwise increase by five-fold of the HIV-specific CTL frequency from 0.04 to 5% increases the probability of detection of HIV-infected cell within 24 h from 0.07 to 0.34 to 0.84 and to 1, respectively. In addition, the model shows that infected motile CD4+ T cells are located faster than non-motile DCs with the probability of detecting them within 24 h increasing from 0.35 to 0.86 and to 1 with a CTL frequency rising from 0.04 to 0.2% and to 1%, respectively.

**Figure 4D** shows the increase of time to locate infected non-motile DCs for a decrease of the T-cell motility from a basal level by 10 and 50% considering an HIV-specific CD8+ T-cell frequency of 1%. If the average T-cell velocity is decreased by 50%, then the probability to locate DCs within 24 h is <0.5. **Figure 4E** depicts a similar dynamics for locating motile infected CD4+ T cells. Note that the motile targets were located within 24 h even with a 50% decrease of the average CD8+ T-cell velocity in all performed simulations.
<figure><img src="image_22.png" /></figure>

DISCUSSION

We have developed a biophysics-based computational model of T-lymphocyte motility that is calibrated using empirical *in vivo* data on T-cell migration in LN tissue representing

three spatial metrics of multicellular systems behavior, i.e., translational speed, turning angle speeds, and meandering index. The model provides the tool to quantify the velocity and the driving force fields in the LN. It enabled us to predict frequency and motility parameters that are required for a
timely detection of productively HIV-infected cells within LNs before they release viral progeny. As such, our study provides a quantitative guide for an informed design of HIV vaccines. Furthermore, as the immunological principles of antigen-specific T-cell activation and immune surveillance imbedded in our model also apply to other infections and cancers, our findings may be used to define the general requirements for any efficient immunotherapeutic intervention against pathogens or cancers in relation to disease-specific parameters and states of lymphoid tissue and T cells. Thus, our model has a significant potential to guide the search for better and more efficient immunotherapies in the near future.

Other processes, e.g., chemotaxis, haptotaxis, and others, can influence the efficacy of target cell search by CTL. The impact of chemotactic migration of T cells toward DCs has been computationally analyzed using a cellular Potts Model (58), an agent-based model (59), and a multicompartmental spatially resolved stochastic model of T-cell circulation (60). The results suggest that the chemoattraction toward target cells modestly speeds up the search process for T cells that successfully find the chemokine-producing DCs. However, a qualitative model presented in (59) suggested that with even weak chemotaxis, substantially lower numbers of CTL are required for sterilizing immunity. Further data-based model-driven research is needed to clarify the contribution of chemotaxis to T-cell migration under normal conditions and during inflammation (61).

Phenomenological Ordinary Differential Equation (ODE) models may also be developed to simulate the interactions between cell populations in the LN. However, these models are not suitable for the present study for three reasons. First, data on T-cell motility in the LN cannot be directly used to calibrate such models, thereby limiting the validity of their predictions. Second, the objective of our study, which is the early detection of HIV-infected T cells and DCs, requires the monitoring of the spatial density of T cells in the LN rather than the total number of T cells. Changes in the spatial distribution of T cells in the LN can be related to spatial mechanisms such as chemotaxis and migration. Therefore, it is crucial to consider spatial aspects in the model. Finally, ODE models based on “mass action”- or “predator-prey”-type parameterizations would require the parameter values specifying a per capita killing rate of target cells. The respective parameter can be determined by the mean time needed for a migrating CTL to locate infected cells. *A priori* estimates of this parameter are not available. It is the spatially resolved model-based simulation that needs to be implemented in order to quantify the killing rate coefficient of the ODE model.

Moving from phenomenological models of spatiotemporal dynamics of immune processes (e.g., the compartmental models, CAMs) to a physics-based description of immune cell migration in complex multicellular tissue environments presents a challenge to mathematical immunology. Advances in the direct visualization of antigen-specific T-cell mobility during their search for and their interaction with antigen-presenting cells within LNs set the basis for diverse modeling approaches (7, 10, 11), which have been so far based on *ad hoc* postulated rules of cell behaviors. Our study gives a biophysics perspective on coordinated cell motility in lymphoid tissues, thus extending the

range of modeling tools available for implementing integrative approaches to the exploration of the immune system.

CPMs have also been applied previously to study intranodal T-cell migration (58). The CPM framework is a valuable tool for a phenomenological description of multicellular patterning, providing realistic simulations of morphological changes for various cell types. The strength of this approach stems from its flexible energetic formalism that allows for extensions to incorporate various biological processes (62). Although the CPM framework has a richer potential for describing individual cell dynamics, including the cell shape, this comes at the expense of (i) a higher-dimensional representation of the cell configuration (e.g., the number of voxels or pixels), (ii) the use of phenomenologically rather than biophysically defined parameters, and (iii) a much higher computational cost to perform simulations required to explore T-cell search strategies. Besides, there is no direct correspondence of most of the CPM parameters with biophysical properties of cells, and the meaning of some CPM parameters is still under debate (62, 63). Moreover, CPM temporal kinetics obtained with the modified Metropolis algorithm does not preserve the detailed balance condition for the underlying stochastic process. This implies that the exact relation between forces of cell interactions and energy terms of CPMs cannot be obtained even for the overdamped dynamics approximation (62, 63).

Computational modeling of multicellular dynamics in lymphoid tissues provides a theoretical tool to be used for a better understanding of the determinants of efficient immune responses against pathogens with a final aim of an optimal manipulation of the immune systems performance (2, 15). Given that the quest for an effective HIV vaccine remains a global priority (64) and that the localization, migration, and frequency of CTLs in LNs determine the extent of virus elimination (17, 19, 36, 56, 65), we sought to use our modeling approach to define threshold frequencies of CTLs in LNs for protection against HIV. Since an HIV infection can influence (i) $CD4^+$ T-cell motility by a direct mechanism involving the HIV Nef protein and (ii) CTL locomotion via an indirect mechanism related to the induction of lymphatic tissue fibrosis, we considered both phenomena to predict the effect of reduction of T-cell motility. We estimated that the frequency of antigen-specific CTL should be about 5% to timely detect and completely eliminate productively infected DCs within 18 h. The time reduces to 4 h for productively infected $CD4^+$ T cells, which are motile. For an HIV-specific T-cell frequency of 1%, we computed that the inhibition of CTL locomotion by two-fold would reduce the probability of detection of infected target cells within 24 h post-infection from 0.84 to 0.42. Thus, the requirements for a prophylactic vaccine for seronegative individuals and an immunotherapeutic intervention of already HIV-infected individuals may differ significantly and are influenced by the state of the lymphatic tissue structure.

Understanding the spatiotemporal dynamics of immune cells globally in the lymphatic system and locally in LNs is considered to be a prerequisite for the development of novel immune interventions in the context of HIV cure strategies (15, 56). To this end, mathematical tools are being increasingly applied to
predict the impacts of trafficking and motility parameters on the efficiency of immune surveillance in health and disease. For example, an optimal surveillance strategy for T cells was analyzed by compartmental modeling of their systemic recirculation and LN transit times using a multicompartmental consideration (66). The protective effect of increased CD4⁺ T-cell trafficking on the dynamics of HIV infection has been recently shown using another compartmental model (67), thus providing a basis for considering cell trafficking as an adjunct therapy option. A multiscale model of *Mycobacterium tuberculosis* infection including an agent-based description of the cellular movement in a two-dimensional simulation grid representing the granuloma was developed and calibrated using non-human primates to derive the prediction of parameters underlying granuloma sterilization (8). However, such modeling attempts are still rare.

In conclusion, the large number of existing mathematical models based on low-resolution descriptions of immune functions has to be further extended and embedded into physiologically distinct compartments and 3D morphological constraints inherent to cells, tissues, and the whole organism. This will then allow the research community not only to get a better quantitative understanding of immune system functioning in infections such as HIV but also enable to build integrative models for antiviral and immunomodulatory drugs of various physical and chemical nature as well as the effects of adoptive cell transfer therapies. We believe that a comprehensive approach to combination therapies based on ART and immunomodulatory drugs affecting a range of processes, including LN fibrosis, the exhaustion of CTLs, and T-cell motility, should rely on formulation and implementation of hybrid spatially resolved multiscale mathematical models of virus infections (8, 9, 68). The here-developed model offers a broadly

applicable generic mathematical tool for linking individual and coordinated cell behaviors that can be used for *in silico* studies to embed immune processes into their spatial context. The physics-based computational model of multicellular dynamics of the immune response in lymphoid tissues provides a solid module that can be universally used in systems immunology studies (2, 6) for the benefit of patients suffering from chronic virus diseases.

## AUTHOR CONTRIBUTIONS

DG, AM, VV, and GB conceived and designed the study. DG, VV, AM, and GB wrote the paper. DG, AB, and NB performed the computational implementation of the model. DG performed the model calibration.

## FUNDING

This work is supported by the Russian Science Foundation (grant 18-11-00171). VV was supported by a grant from the RUDN University Program 5-100. AM is also supported by a grant from the Spanish Ministry of Economy, Industry and Competitiveness and FEDER grant no. SAF2016-75505-R (AEI/MINEICO/FEDER, UE) and the María de Maeztu Programme for Units of Excellence in R&D (MDM-2014-0370).

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2019.01213/full#supplementary-material

## REFERENCES

1. Munoz MA, Biro M, Weninger W. T cell migration in intact lymph nodes in vivo. *Curr Opin Cell Biol.* (2014) 30:17-24. doi: 10.1016/j.ceb.2014.05.002

2. Wong HS, Germain RN. Robust control of the adaptive immune system. *Semin Immunol.* (2017) 36:17-27. doi: 10.1016/j.smin.2017.12.009

3. Ludewig B, Stein JV, Sharpe J, Cervantes-Barragan L, Thiel V, Bocharov G. A global "imaging" view on systems approaches in immunology. *Eur J Immunol.* (2012) 42:3116-25. doi: 10.1002/eji.201242508

4. Worbs T, Förster R. T cell migration dynamics within lymph nodes during steady state: an overview of extracellular and intracellular factors influencing the basal intranodal T cell motility. In: Dustin M, McGovern D, editors. *Visualizing Immunity Current Topics in Microbiology and Immunology*. Berlin: Heidelberg: Springer (2009). p. 71-105. doi: 10.1007/978-3-540-93864-4_4

5. Mandl JN, Liou R, Klauschen F, Vrisekoop N, Monteiro JP, Yates AJ, et al. Quantification of lymph node transit times reveals differences in antigen surveillance strategies of naïve CD4⁺ and CD8⁺ T cells. *Proc Natl Acad Sci USA.* (2012) 109:18036-41. doi: 10.1073/pnas.1211717109

6. Germain RN. Will systems biology deliver its promise and contribute to the development of new or improved vaccines? What really constitutes the study of "systems biology" and how might such an approach facilitate vaccine design. *Cold Spring Harb Perspect Biol.* (2018) 10:a033308. doi: 10.1101/cshperspect.a033308

7. Gong C, Mattila JT, Miller M, Flynn JL, Linderman JJ, Kirschner D. Predicting lymph node output efficiency using systems biology. *J Theor Biol.* (2013) 335:169-84. doi: 10.1016/j.jtbi.2013.06.016

8. Pienaar E, Matern WM, Linderman JJ, Bader JS, Kirschner DE. Multi-scale model of *Mycobacterium tuberculosis* infection maps metabolite and gene perturbations to granuloma sterilization predictions. *Infect Immun.* (2016) 84:1650-69. doi: 10.1128/IAI.01438-15

9. Bouchnita A, Bocharov G, Meyerhans A, Volpert V. Hybrid approach to model the spatial regulation of T cell responses. *BMC Immunol.* (2017) 18:29. doi: 10.1186/s12865-017-0205-0

10. Beauchemin C, Dixit NM, Perelson AS. Characterizing T cell movement within lymph nodes in the absence of antigen. *J Immunol.* (2007) 178:5505-12. doi: 10.4049/jimmunol.178.9.5505

11. Beltman JB, Marée AFM, Lynch JN, Miller MJ. Boer RJ de Lymph node topology dictates T cell migration behavior. *J Exp Med.* (2007) 204:771-80. doi: 10.1084/jem.20061278

12. Drasdo D, Hoehme S, Block M. On the role of physics in the growth and pattern formation of multi-cellular systems: what can we learn from individual-cell based models? *J Stat Phys.* (2007) 128:287. doi: 10.1007/s10955-007-9289-x

13. Risler T. Cytoskeleton and cell motility. In: Meyers RA, editor. *Encyclopedia of Complexity and Systems Science*. New York, NY: Springer (2013). p. 1-53. doi: 10.1007/978-3-642-27737-5_112-3

14. Camley BA, Rappel W-J. Physical models of collective cell motility: from cell to tissue. *J Phys Appl Phys.* (2017) 50:113002. doi: 10.1088/1361-6463/aa56fe

15. Krummel MF, Bartumeus F, Gérard A. T cell migration, search strategies and mechanisms. *Nat Rev Immunol.* (2016) 16:193-201. doi: 10.1038/nri.2015.16
16. Gadhamsetty S, Marée AFM, Boer DJR, Beltman JB. Tissue dimensionality influences the functional response of cytotoxic T lymphocyte-mediated killing of targets. *Front Immunol*. (2017) 7:668. doi: 10.3389/fimmu.2016.00668

17. Foley MH, Forcier T, McAndrew E, Gonzalez M, Chen H, Juelg B, et al. High avidity CD8⁺ T cells efficiently eliminate motile HIV-infected targets and execute a locally focused program of anti-viral function. *PLoS ONE*. (2014) 9:e87873. doi: 10.1371/journal.pone.0087873

18. Murooka TT, Deruaz M, Marangoni F, Vrbanac VD, Seung E, von Andrian UH, et al. HIV-infected T cells are migratory vehicles for viral dissemination. *Nature*. (2012) 490:283–7. doi: 10.1038/nature11398

19. Fukazawa Y, Park H, Cameron MJ, Lefebvre F, Lum R, Coombes N, et al. Lymph node T cell responses predict the efficacy of live attenuated SIV vaccines. *Nat Med*. (2012) 18:1673–81. doi: 10.1038/nm.2934

20. Mohammadi P, Desfarges S, Bartha I, Joos B, Zangger N, Muñoz M, et al. 24 hours in the life of HIV-1 in a T cell line. *PLoS Pathog*. (2013) 9:e1003161. doi: 10.1371/journal.ppat.1003161

21. O'Neill ME. PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation. Claremont, CA: Harvey Mudd College (2014).

22. Read MN, Bailey J, Timmis J, Chtanova T. Leukocyte motility models assessed through simulation and multi-objective optimization-based model selection. *PLoS Comput Biol*. (2016) 12:e1005082. doi: 10.1371/journal.pcbi.1005082

23. Zimmermann J, Camley BA, Rappel W-J, Levine H. Contact inhibition of locomotion determines cell–cell and cell–substrate forces in tissues. *Proc Natl Acad Sci USA*. (2016) 113:2660–5. doi: 10.1073/pnas.152230113

24. Kitano M, Yamazaki C, Takumi A, Ikeno T, Hemmi H, Takahashi N, et al. Imaging of the cross-presenting dendritic cell subsets in the skin-draining lymph node. *Proc Natl Acad Sci USA*. (2016) 113:1044–9. doi: 10.1073/pnas.1513607113

25. Hakim V, Silberzan P. Collective cell migration: a physics perspective. *Rep Prog Phys*. (2017) 80:076601. doi: 10.1088/1361-6633/aa65ef

26. Basan M, Prost J, Joanny J-F, Elgeti J. Dissipative particle dynamics simulations for biological tissues: rheology and competition. *Phys Biol*. (2011) 8:026014. doi: 10.1088/1478-3975/8/2/026014

27. Liu MB, Liu GR, Zhou LW, Chang JZ. Dissipative particle dynamics (DPD): an overview and recent developments. *Arch Comput Methods Eng*. (2015) 22:529–56. doi: 10.1007/s11831-014-9124-x

28. Bell GI. Models for the specific adhesion of cells to cells. *Science*. (1978) 200:618–27. doi: 10.1126/science.347575

29. Kong F, García AJ, Mould AP, Humphries MJ, Zhu C. Demonstration of catch bonds between an integrin and its ligand. *J Cell Biol*. (2009) 185:1275–84. doi: 10.1083/jcb.200810002

30. Basu R, Huse M. Mechanical communication at the immunological synapse. *Trends Cell Biol*. (2017) 27:241–54. doi: 10.1016/j.tcb.2016.10.005

31. Segel GB, Cokelet GR, Lichtman MA. The measurement of lymphocyte volume: importance of reference particle deformability and counting solution tonicity. *Blood*. (1981) 57:894–9.

32. Wiśniewski JR, Hein MY, Cox J, Mann M. A “proteomic ruler” for protein copy number and concentration estimation without spike-in standards. *Mol Cell Proteomics*. (2014) 13:3497–506. doi: 10.1074/mcp.M113.037309

33. Illmer P, Erlebach C, Schinner F. A practicable and accurate method to differentiate between intra- and extracellular water of microbial cells. *FEMS Microbiol Lett*. (1999) 178:135–9. doi: 10.1111/j.1574-6968.1999.tb13769.x

34. Brown GC. Total cell protein concentration as an evolutionary constraint on the metabolic control distribution in cells. *J Theor Biol*. (1991) 153:195–203. doi: 10.1016/S0022-5193(05)80422-9

35. Turgeon ML. *Clinical Hematology: Theory and Procedures*. Philadelphia, PA: Lippincott Williams & Wilkins (2005).

36. Lim TS, Goh JKH, Mortellaro A, Lim CT, Hämmerling GJ, Ricciardi-Castagnoli P. CD80 and CD86 differentially regulate mechanical interactions of T-cells with antigen-presenting dendritic cells and B-cells. *PLoS ONE*. (2012) 7:e45185. doi: 10.1371/journal.pone.0045185

37. Miller MJ, Wei SH, Cahalan MD, Parker I. Autonomous T cell trafficking examined *in vivo* with intravital two-photon microscopy. *Proc Natl Acad Sci USA*. (2003) 100:2604–9. doi: 10.1073/pnas.2628040100

38. Zhang ML, Guo AX, VandenBussche CJ. Morphologists overestimate the nuclear-to-cytoplasmic ratio. *Cancer Cytopathol*. (2016) 124:669–77. doi: 10.1002/cncy.21735

39. Basu R, Whitlock BM, Husson J, Le Floc'h A, Jin W, Oyler-Yaniv A, et al. Cytotoxic T cells use mechanical force to potentiate target cell killing. *Cell*. (2016) 165:100–10. doi: 10.1016/j.cell.2016.01.021

40. Purcell EM. Life at low Reynolds number. *Am J Phys*. (1977) 45:3–11. doi: 10.1119/1.10903

41. Friedl P, Bröcker E-B. The biology of cell locomotion within three-dimensional extracellular matrix. *Cell Mol Life Sci*. (2000) 57:41–64. doi: 10.1007/s000180050498

42. Hons M, Kopf A, Hauschild R, Leithner A, Gaertner F, Abe J, et al. Chemokines and integrins independently tune actin flow and substrate friction during intranodal migration of T cells. *Nat Immunol*. (2018) 19:606–16. doi: 10.1038/s41590-018-0109-z

43. Goya GF, Marcos-Campos I, Fernández-Pacheco R, Sáez B, Godino J, Asín L, et al. Dendritic cell uptake of iron-based magnetic nanoparticles. *Cell Biol Int*. (2008) 32:1001–5. doi: 10.1016/j.cellbi.2008.04.001

44. Han X, Kloeden PE. *Random Ordinary Differential Equations and Their Numerical Solution*. Singapore: Springer (2017). doi: 10.1007/978-981-10-6265-0

45. Rothoeft T, Balkow S, Krummen M, Beissert S, Varga G, Loser K, et al. Structure and duration of contact between dendritic cells and T cells are controlled by T cell activation state. *Eur J Immunol*. (2006) 36:3105–17. doi: 10.1002/eji.200636145

46. Mempel TR, Henrickson SE. Andrian UH von. T-cell priming by dendritic cells in lymph nodes occurs in three distinct phases. *Nature*. (2004) 427:154. doi: 10.1038/nature02238

47. Miller MJ, Hejazi AS, Wei SH, Cahalan MD, Parker I. T cell repertoire scanning is promoted by dynamic dendritic cell behavior and random T cell motility in the lymph node. *Proc Natl Acad Sci USA*. (2004) 101:998–1003. doi: 10.1073/pnas.0306407101

48. Blattman JN, Antia R, Sourdive DJD, Wang X, Kaech SM, Murali-Krishna K, et al. Estimating the precursor frequency of naïve antigen-specific CD8 T Cells. *J Exp Med*. (2002) 195:657–64. doi: 10.1084/jem.20001021

49. Bastidas S, Graw F, Smith MZ, Kuster H, Günthard HF, Oxenius A. CD8⁺ T cells are activated in an antigen-independent manner in HIV-infected individuals. *J Immunol*. (2014) 192:1732–44. doi: 10.4049/jimmunol.1302027

50. Haynes BF, Shaw GM, Korber B, Kelsoe G, Sodroski J, Hahn BH, et al. HIV-host interactions: implications for vaccine design. *Cell Host Microbe*. (2016) 19:292–303. doi: 10.1016/j.chom.2016.02.002

51. Grossman Z, Feinberg MB, Paul WE. Multiple modes of cellular activation and virus transmission in HIV infection: a role for chronically and latently infected cells in sustaining viral replication. *Proc Natl Acad Sci USA*. (1998) 95:6314–9. doi: 10.1073/pnas.95.11.6314

52. Grossman Z, Meier-Schellersheim M, Paul WE, Picker LJ. Pathogenesis of HIV infection: what the virus spares is as important as what it destroys. *Nat Med*. (2006) 12:289–95. doi: 10.1038/nm1380

53. Cheynier R, Henrichwark S, Hadida F, Pelletier E, Oksenhendler E, Autran B, et al. HIV and T cell expansion in splenic white pulps is accompanied by infiltration of HIV-specific cytotoxic T lymphocytes. *Cell*. (1994) 78:373–87. doi: 10.1016/0092-8674(94)90417-0

54. Haase AT. Population biology of HIV-1 infection: viral and CD4⁺ T cell demographics and dynamics in lymphatic tissues. *Annu Rev Immunol*. (1999) 17:625–56. doi: 10.1146/annurev.immunol.17.1.625

55. Kityo C, Makamdp KN, Rothenberger M, Chipman JG, Hoskulsson T, Beilman GJ, et al. Lymphoid tissue fibrosis is associated with impaired vaccine responses. *J Clin Invest*. (2018) 128:2763–73. doi: 10.1172/JCI97377

56. Dimopoulos Y, Moysi E, Petrovas C. The lymph node in HIV pathogenesis. *Curr HIV/AIDS Rep*. (2017) 14:133–40. doi: 10.1007/s11904-017-0359-7

57. Fackler OT, Murooka TT, Imle A, Mempel TR. Adding new dimensions: towards an integrative understanding of HIV-1 spread. *Nat Rev Microbiol*. (2014) 12:563–74. doi: 10.1038/nrmicro3309

58. Vroomans RMA, Marée AFM, de Boer RJ, Beltman JB. Chemotactic migration of T cells towards dendritic cells promotes the detection of rare antigens. *PLoS Comput Biol*. (2012) 8:e1002763. doi: 10.1371/journal.pcbi.1002763

59. Kadolsky UD, Yates AJ. How is the effectiveness of immune surveillance impacted by the spatial distribution of spreading infections? *Philos Trans R Soc B Biol Sci*. (2015) 370:20140289. doi: 10.1098/rstb.2014.0289
60. Textor J, Henrickson SE, Mandl JN, von Andrian UH, Westermann J, de Boer RJ, et al. Random migration and signal integration promote rapid and robust T cell recruitment. *PLoS Comput Biol.* (2014) 10:e1003752. doi: 10.1371/journal.pcbi.1003752

61. Banigan EJ, Harris TH, Christian DA, Hunter CA, Liu AJ. Heterogeneous CD8⁺ T cell migration in the lymph node in the absence of inflammation revealed by quantitative migration analysis. *PLoS Comput Biol.* (2015) 11:e1004058. doi: 10.1371/journal.pcbi.1004058

62. Scianna M, Preziosi L. Cellular Potts Models: Multiscale Extensions and Biological Applications. New York, NY: CRC Press (2013). doi: 10.1201/b14075

63. Anderson ARA, Chaplain MAJ, Rejniak KA (eds.). *Single-Cell-Based Models in Biology and Medicine*. Basel: Birkhäuser Basel (2007). doi: 10.1007/978-3-7643-8123-3

64. Haynes BF, Burton DR. Developing an HIV vaccine. *Science*. (2017) 355:1129-30. doi: 10.1126/science.aan0662

65. Connick E, Folkvord JM, Lind KT, Rakasz EG, Miles B, Wilson NA, et al. Compartmentalization of simian immunodeficiency virus replication within secondary lymphoid tissues of rhesus macaques is linked to disease stage and inversely related to localization of virus-specific CTL. *J Immunol.* (2014) 193:5613-25. doi: 10.4049/jimmunol.1401161

66. Lee M, Mandl JN, Germain RN, Yates AJ. The race for the prize: T cell trafficking strategies for optimal surveillance. Blood. (2012) 120:1432-8. doi: 10.1182/blood-2012-04-424655

67. Fryer HR, Wolinsky SM, McLean AR. Increased T cell trafficking as adjunct therapy for HIV-1. *PLoS Comput. Biol.* (2018) 14:e1006028. doi: 10.1371/journal.pcbi.1006028

68. Bocharov G, Volpert V, Ludewig B, Meyerhans A. *Mathematical Immunology of Virus Infections*. Cham: Springer (2018). doi: 10.1007/978-3-319-72317-4

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Grebennikov, Bouchnita, Volpert, Bessonov, Meyerhans and Bocharov. 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.