Simulation Experiments in Daphne¶

Modeling the inter-zonal dynamics of germinal center response¶

Previous work has experimentally demonstrated bidirectional germinal center (GC) B cell movement between zones [Victora2010], [Allen2007], which supports the iterative cycles of mutation and selection (Kepler and Perelson, 1993b, Kepler and Perelson, 1993a, [Meyer_Hermann2001] ) necessary for effective affinity maturation. In this work, we try to 1) generate a computational model to further understand the details of the interzonal migration of germinal center B cells following the experimental work done by Victora et al; 2) test and improve our agent-based modeling software, Daphne, by running simulations of the above-mentioned inter-zonal dynamics model. Ordinary Differential Equation model of inter-zonal dynamics

Ordinary Differential Equation model of inter-zonal dynamics¶

We first develop the ordinary different equation (ODE) model following the work done by Victora, et al (Figure 1 and Figure 2).

Figure 1: Interzonal migration in DCs. (A) experimental protocol. (B) Time series showing position of individual cells photoactivated in the DZ. (C) Time series for cells photoactivated in LZ. (Figure was taken from [Victora2010])

Figure 2: Quantification summary of seven independent experiments from Victora, et al.

Based on the experiments in Victora’s work as well as the current knowledge about the GC response, we set up the following model (Figure 3).

Figure 3: Interzonal dynamics model.

The following ODEs were developed for the model,

(1)\[\begin{split} \frac{dB}{dt} &= \rho B - \mu B + \nu C \\ \frac{dC}{dt} &= \mu B - \nu C - \gamma C - \alpha C \\ \frac{dA}{dt} &= \alpha C - \beta A\end{split}\]

Daphne model of inter-zonal dynamics¶

Based on the ODE model built in the above sections, a model was set up in Daphne. The model contains the following components,

Germinal center space. This was a cubic space with a dimension of 100 μm x 100 μm x 100 μm (Figure 4).
Dark zone. The space was contained in the germinal center cube and was defined by the cytokine CXCL12, which followed a Gaussian distribution with a center at (75, 50, 50) (Figure 4).
Light zone. The space was contained in the germinal center cube and was defined by the cytokine CXCL13, which followed a Gaussian distribution with a center at (25, 50, 50) (Figure 4).

Modeling centroblast and centrocyte populations¶

Centroblast. The B cell population expresses cell surface receptor CXCR4 that can bind to CXCL12. Formation of CXCL12 and CXCR4 complex leads to activation of chemotaxis that drives migration following the gradient as well as cell division and differentiation.
Centrocyte. The B cell population expresses cell surface receptor CXCR5 that can bind to CXCL13. Formation of CXCL13 and CXCR5 complex leads to activation of chemotaxis that drives migration following the gradient as well as cell death and differentiation.

This simplified model did not include receptor dynamics - internalization and degradation of receptors.

Figure 4: Germinal center model in Daphne.

Results¶

ODE model fitting¶

The experimental data were extracted from the work done by Victora, et al. (Fig. 1). Then the ODE model was fit to data in order to estimate the parameter listed in Fig 3. A Markov Chain Monte Carlo method with the Metropolis-Hasting algorithm was run for estimation. The results are listed below,

\(\rho\) =0.0319, \(\gamma\) =0.0262
\(\mu\) =0.1175, \(\nu\) =0.923469
\(\alpha\) =0.318, \(\beta\) =0.06502
\(B_0\) =50, \(C_0\) =7.12, \(A_0\) =66.7

The fitted values are plotted in Fig 5.

Figure 5: Model fitting of ODE model. Dashed lines are experiment data and solid lines are fitted values.

Simulation in Daphne¶

With the parameters estimated using the ODE model, we built the Daphne model. Two experiments were carried out, DZ to LZ and LZ to DZ migration. The simulation was run 6 hours and B cells movement in the DZ or LZ were recorded, as in the original reference paper. The (stochastic) simulation was repeated at least 5 times. Results are shown in Fig 6.

Figure 6: Simulation of inter-zonal dynamics in Daphne. The blue and green lines are Daphne simulation results. The red dashed lines are experimental data and green solid lines are ODE results.

Lymphocyte locomotion kinetics: Comparison of 2D ex vivo and 3D intravital preparations.¶

The aim of the present study was to determine whether or not the results of the 2D ex vivo lymphocyte studies carried out by the Haugh group are likely to provide information relevant to understanding lymphocyte motion in vivo.

To address this question, we collected data from murine T cells in the Haugh 2D system and from murine GC T cells in the Haberman 2p intravital system. These data were generated prior to our undertaking this comparison; they were not generated specifically to answer the current question. We used T cells in this study because the T cell system is more mature than the B cell system in the Haugh laboratory and the data are more reliable for that reason. We intend to repeat the analysis for B cells.

We analyzed the data by fitting them to a first-order Langevin process of appropriate spatial dimension. The parameters \(\gamma\) and \(\zeta\) were fit to each trajectory individually, as were the error variances and initial values for the positions and velocities.

The Langevin model is given by the Ito differential equation

(2)\[\begin{split} dX &= V dt \\ dV &= -\gamma V dt + \zeta \sqrt{\gamma dW}\end{split}\]

where \(X\) is the position of the cell, and \(V\) is the cell velocity. \(\gamma\) is the inverse of the persistence time, and \(\zeta\) is proportional to the mean speed (the constant of proportionality depends on the spatial dimension).

The model is fit by deriving the covariance function and mean functions from Eq., and maximizing the likelihood over \(\gamma, \zeta, X(0), V(0)\) and \(\sigma\), where \(\sigma\) is the measurement error.

The instantaneous speeds and direction vectors are estimated directly from the most-probable trajectory under the posterior distribution.

The sampling interval for the 3D intravital data was 25 seconds, while that for the 2D ex vivo data was 3 seconds. For this analysis we undersampled each 2D time series by six-fold to achieve a sampling interval close to that used for the 3D data, while maintaining sufficiently long trajectories.

We required 60 points per trajectory for analysis, and accepted all trajectories from both data sets that met that criterion. The 2D data had 7 trajectories (out of 10) that met the criterion; the 3D data had 48 trajectories (out of over 2000) that met the criterion.

Statistics¶

For each trajectory, we estimated the instantaneous speed and direction vector at each time in the usual way

(3)\[ s(t) = \sqrt{\sum_{i=1}^{\Delta} V_i^2(t)}\]

where \(\Delta\) is the dimension of the relevant space (2 or 3), and

(4)\[ u(t) = s^{-1}V(t)\]

is the unit direction vector.

The instantaneous reorientation rate is defined as the magnitude of the rate of change of the direction vector

(5)\[ \rho(t) = \| \frac{d u(t)}{dt} \|\]

and is estimated by the finite difference

(6)\[ \rho(t) = \| \frac{u(t_{i+1}) - u(t_i)}{t_{i+1}-t_i} \|\]

The speed autocorrelation function is estimated as

(7)\[ f_s(\tau) = \frac{\sum s(t_i) s(t_i + \tau)}{\sum s(t_i)^2}\]

and the direction autocorrelation function is

(8)\[ f_s(\tau) = \frac{\sum u(t_i) \cdot u(t_i + \tau)}{\sum \| u(t_i) \|^2}\]

In both cases, the parameter \(\tau\) is known as the lag.

Results¶

Figure 7 shows the fitted trajectories for the largest three single-cell data sets from each group. In both cases, it appears that the trajectories turn during periods of low speed. The plots are presented for reference only. Although the speeds and the scale at which speeds and directions change appear to be similar, more direct comparison is required for verification.

Figure 7: Trajectories from the 3D intravital (top row) and 2D ex vivo systems. The three longest trajectories for each type are shown. The speed is indicated by color going from blue (slow) to red (fast)..

The overall mean speeds (root mean-square speeds) for the two datasets are remarkably similar given the enormous differences between their treatments and environments, with the mean speed of the intravital trajectories about 7% smaller than that of the ex vivo trajectories:

3D intravital: 5.28 m/min 2D ex vivo: 5.70 m/min

To examine the turning kinetics more closely, we plot the reorientation rate against the speed for the same three longest trajectories in either group. Figure 8 shows these scatter plots, which show that the scatter patterns of trajectories in both sets are very similar. Computation of the correlation coefficient between the log reorientation rate and the log speed for each trajectory yields the following table:

Figure 8: The reorientation rate vs the speed for each of the six trajectories shown in Figure 7.

Although the intravital trajectories show a somewhat stronger anticorrelation between reorientation rate and speed than do the ex vivo trajectories, the most striking result is the consistency of the strong anticorrelation in each case.

The estimated Langevin parameters are themselves informative about the characteristics of the trajectories. Figure 9. shows scatter plots of the estimated values of \(\gamma\) and \(\zeta\) for all of the trajectories analyzed. The two scatter plots are completely overlapping. The 95% percent confidence ellipses also overlap, though the inclusion of more data in the ex vivo data sets might render them more distinct.

Figure 9: Estimated zeta vs estimated gamma (filled circles) for all of the 2D (red) and 3D (blue) trajectories. The diamonds with show the position and 95% confidence intervals for the means of gamma and zeta.

Finally, we examine the autocorrelation functions for speed and direction. The widths of these functions are determined by the rate at which the cell loses memory of its speed or direction, respectively. These functions are key reporters of the cell’s local behavior.

Figure 10, bottom row, shows comparisons between the autocorrelation functions for the intravital and ex vivo trajectories. It is remarkable how similar the initial declines in autocorrelation are between the two sets. Beyond about half a minute, however the intravital correlation declines at a significantly greater rate. We hypothesize that the observed difference in autocorrelation is due to cellular crowding in the intravital system. In the germinal center, B cells are very tightly packed. It is believed that the proportion of fluorescent T cells among all lymphocytes in the field of vision in this intravital preparation is small. In contrast, the ex vivo preps are isolated single cells.

Figure 10: The autocorrelation functions for speed (left) and direction vectors (right). The top row refers to simulated data. Sparse and crowded are defined in the text. The lower row refers to experimental data as indicated by color.

To investigate this hypothesis, we simulated the effect of cellular crowding using the Plaza Sur simulation platform. We used toroidal boundary conditions: a cell passing through any boundary of the simulation cube reenters at the corresponding point on the opposite face with identical velocity. We assumed mechanical repulsion among cells in contact with each other. Cells closer than required for surface-surface contact experience a force proportional to the penetration distance. We omitted any adhesive forces among cells.

For each run, we used 500 cells, and varied the total size of the tissue volume from 109 to 1.25 × 105 m3. At the former volume, cells are very unlikely to encounter each other. At the latter volume, cells are squeezed down to one half of their own natural volume. For each volume, we computed the speed and direction autocorrelations. Figure 11 shows the variation of the direction autocorrelation as the cells become more crowded. For cells packed more densely than 1 cell per 20m cube (labeled “5” in the plot) there is a steady steepening of the autocorrelation function.

Figure 11: Direction autocorrelation functions for simulated cell populations of decreasing total volume. The total number of cells is 500 in each run, and the volume decreases by a factor of two between adjacent curves.

The upper row of Figure 10 shows a comparison of the autocorrelation functions for a very sparse system compared to those for a dense system. The packing fraction is the total cellular volume divided by the volume of the tissue containing them.

Sparse: 500 cells in a tissue volume of 1 × 109 m3 for a packing fraction of 2.5× 10-4.
Crowded: 500 cells in a tissue volume of 2.5 × 105 m3 for a packing fraction of about 1.0.

Discussion¶

We have analyzed data from the Haugh group’s 2D ex vivo system, and data from the Haberman group’s intravital 2p system and compared the kinetics of cellular motion in the two cases. The similarities between the cell kinetics in the two cases are striking. In spite of the differences in geometry and environment, the locomotive behaviors are remarkable in their sameness.

Where the two systems do have distinct behaviors—in the autocorrelations—our simulations suggest that these differences can be accounted for by cell-cell mechanical contact, which occur at very high frequency in the in the intravital system but are purposely excluded from the ex vivo system.

We believe these findings justify the continued use of the ex vivo system to probe lymphocyte locomotion and to extend these studies to morphodynamics and chemotaxis.