Document Detail

 A Model for Migratory B Cell Oscillations from Receptor Down-Regulation Induced by External Chemokine Fields. Jump to Full Text MedLine Citation: PMID:  23296998     Owner:  NLM     Status:  Publisher Abstract/OtherAbstract: A long-standing paradigm in B cell immunology is that effective somatic hypermutation and affinity maturation require cycling between the dark zone and light zone of the germinal center. The cyclic re-entry hypothesis was first proposed based on considerations of the efficiency of affinity maturation using an ordinary differential equations model for B cell population dynamics. More recently, two-photon microscopy studies of B cell motility within lymph nodes in situ have revealed the complex migration patterns of B lymphocytes both in the preactivation follicle and post-activation germinal center. There is strong evidence that chemokines secreted by stromal cells and the regulation of cognate G-protein coupled receptors by these chemokines are necessary for the observed spatial cell distributions. For example, the distribution of B cells within the light and dark zones of the germinal center appears to be determined by the reciprocal interaction between the level of the CXCR4 and CXCR5 receptors and the spatial distribution of their respective chemokines CXCL12 and CXCL13. Computer simulations of individual-based models have been used to study the complex biophysical and mechanistic processes at the individual cell level, but such simulations can be challenging to parameterize and analyze. In contrast, ordinary differential equations are more tractable, but traditional compartment model formalizations ignore the spatial chemokine distribution that drives B cell redistribution. Motivated by the desire to understand the motility patterns observed in an individual-based simulation of B cell migration in the lymph node, we propose and analyze the dynamics of an ordinary differential equation model incorporating explicit chemokine spatial distributions. While there is experimental evidence that B cell migration patterns in the germinal center are driven by extrinsically regulated differentiation programs, the model shows, perhaps surprisingly, that feedback from receptor down-regulation induced by external chemokine fields can give rise to spontaneous interzonal and intrazonal oscillations in the absence of any extrinsic regulation. While the extent to which such simple feedback mechanisms contributes to B cell migration patterns in the germinal center is unknown, the model provides an alternative hypothesis for how complex B cell migration patterns might arise from very simple mechanisms. Authors: Cliburn Chan; Matthew Billard; Samuel A Ramirez; Harald Schmidl; Eric Monson; Thomas B Kepler Related Documents : 22951908 - Glutathione peroxidase 1 activity dictates the sensitivity of glioblastoma cells to oxi...22931348 - Effects of a cell-imprinted poly(dimethylsiloxane) surface on the cellular activities o...22992018 - Metabolism of inorganic arsenic in intestinal epithelial cell lines.23722198 - Induction of apoptosis by high-dose gold nanoparticles in nasopharyngeal carcinoma cells.17065808 - Activation of kinases upon volume changes: role in cellular homeostasis.3383578 - Quantitative effects of short- and long-term vasectomy on mouse spermatogenesis and spe... Publication Detail: Type:  JOURNAL ARTICLE     Date:  2013-1-8 Journal Detail: Title:  Bulletin of mathematical biology     Volume:  -     ISSN:  1522-9602     ISO Abbreviation:  Bull. Math. Biol.     Publication Date:  2013 Jan Date Detail: Created Date:  2013-1-8     Completed Date:  -     Revised Date:  - Medline Journal Info: Nlm Unique ID:  0401404     Medline TA:  Bull Math Biol     Country:  - Other Details: Languages:  ENG     Pagination:  -     Citation Subset:  - Affiliation: Department of Biostatistics and Bioinformatics, Duke University Medical Center, Durham, NC, 27705, USA, cliburn.chan@duke.edu. Export Citation: APA/MLA Format     Download EndNote     Download BibTex MeSH Terms Descriptor/Qualifier:

From MEDLINE®/PubMed®, a database of the U.S. National Library of Medicine

 Full Text Journal Information Journal ID (nlm-ta): Bull Math Biol Journal ID (iso-abbrev): Bull. Math. Biol ISSN: 0092-8240 ISSN: 1522-9602 Publisher: Springer-Verlag, New York Article Information Download PDF © The Author(s) 2012 Received Day: 14 Month: 5 Year: 2012 Accepted Day: 15 Month: 11 Year: 2012 Electronic publication date: Day: 8 Month: 1 Year: 2013 pmc-release publication date: Day: 8 Month: 1 Year: 2013 Print publication date: Month: 1 Year: 2013 Volume: 75 Issue: 1 First Page: 185 Last Page: 205 PubMed Id: 23296998 ID: 3547247 Publisher Id: 9799 DOI: 10.1007/s11538-012-9799-9
 A Model for Migratory B Cell Oscillations from Receptor Down-Regulation Induced by External Chemokine Fields Cliburn ChanAff1 Address: cliburn.chan@duke.edu Matthew BillardAff2 Samuel A. RamirezAff3 Harald SchmidlAff1 Eric MonsonAff4 Thomas B. KeplerAff5 Department of Biostatistics and Bioinformatics, Duke University Medical Center, Durham, NC 27705 USA Thurston Arthritis Research Center, The University of North Carolina at Chapel Hill, Chapel Hill, NC 27599 USA Program in Computational Biology and Bioinformatics, Duke University, Durham, NC 27710 USA Duke University Visualization Technology Group, Duke University, Durham, NC 27708 USA Department of Microbiology, Boston University School of Medicine, 72E Concord St, Boston, MA 02118 USA

Introduction

The evolution of high-affinity specific antibodies by long-lived B cells is driven by a process known as affinity maturation that occurs in the germinal center of lymph nodes. In this process, the germinal center (GC) is partitioned into a dark zone (DZ), consisting largely of rapidly dividing B cells known as centroblasts, and a light zone (LZ), consisting largely of B cells known as centrocytes interacting with follicular dendritic cells (FDC). It is believed that somatic hypermutation which introduces random changes in the antibody nucleotide sequence occurs within centroblasts in the DZ, while centrocytes in the LZ interact with and compete for immune complexes bound to FDC (Allen et al. 2007a; Shlomchik and Weisel 2012). A long-held hypothesis of cyclic re-entry is that the periodic migration of B cells from the DZ to the LZ and vice versa is critical for the efficiency of affinity maturation (Kepler and Perelson 1993; Kepler et al. 1993; Meyer-Hermann et al. 2001). FDCs present antigen bound on Fc receptor-captured antibodies on their cell surface, and centrocytes compete for binding to these antigens. Centrocytes with high affinity B cell receptors are more likely to successfully bind antigen and receive survival signals, while centrocytes with low affinity receptors fail to bind and undergo apoptosis. Successful centrocytes may then reenter the DZ for proliferation and another iteration of selection, or exit the germinal center as memory B cells or long-lived plasma cells.

How B cells migrate in the lymph node is hence critical for understanding the generation of high affinity long-lived memory and plasma cells that are the basis of humoral immunity. Naive B cells are believed to be attracted to the preactivation follicle primarily by the chemokine CXCL13, although lipid ligands that bind to the EBI2 receptor and CCR7:CCL19/CCL21 receptor-ligand interactions also modulate the B cell spatial distribution in the follicle (Gatto et al. 2011). In the post-activation germinal center, the migration of B cells between the DZ and LZ is driven by the chemokines CXCL12 found mainly in the LZ and CXCL13 in the DZ. These chemokines are recognized by the G-protein coupled receptors (GPCR) CXCR4 and CXCR5, with CXCR4 binding to CXCL12 and CXCR5 binding to CXCL13 (Allen et al. 2004, 2007a).

Pioneering work by Sally Zigmond has described receptor internalization (and resulting loss of sensitivity to a chemokine gradient) as an important aspect of GPCR-mediated chemotaxis (Zigmond 1981; Zigmond et al. 1982). Estimated receptor levels are in the 104 range (10,000–50,000) of receptors per cell (Zigmond 1981). Upon ligand binding, GPCRs signal to G-protein, become phosphorylated by GPCR kinase (effectively desensitizing the receptor by dissociating G-protein subunits), and internalize via one of two major pathways. One internalization pathway is fast and involves clathrin-mediated endocytosis. The other is slower and uses a lipid-raft/caveolae pathway. The clathrin pathway involves the recruitment of arrestin to the receptor, which can act as a scaffold for further signaling events. Receptors internalized in either pathway can potentially be recycled, or degraded. How chemokine receptors respond to the local chemokine field over time is hence likely to be a major regulatory mechanism for the migration behavior of B cells. Indeed, the literature describes alterations in chemokine receptor expression balance as the fundamental basis for directional migration within the lymph node and germinal center (Allen et al. 2004; Hardtke et al. 2005; Reif et al. 2002).

With the advent of two-photon microscopy, we can now observe individual B cell dynamics in situ within a developing germinal center (Allen et al. 2007b; Schwickert et al. 2007; Hauser et al. 2007b; Victora et al. 2010). However, two-photon microscopy is restricted to the visualization of relatively small regions and short time-spans. Computational modeling is therefore a valuable adjunct for inference beyond these short time and space scales, providing mechanistic insight into long range/long duration phenomena such as the relationship between B cell migration patterns and the efficacy of somatic hypermutation (Kleinstein 2002; Meyer-Hermann et al. 2009; Figge and Meyer-Hermann 2011). As traditional ordinary differential equation (ODE) models ignore the spatial inhomogeneity of the chemokine fields, computational simulations of individual-based models (IBM) may be more appropriate vehicles for understanding how emergent behavior arises from the interactions of single B cells with their environment and other cells (Figge 2005; Bogle and Dunbar 2009; Germain et al. 2011; Beltman et al. 2011).

The detailed biology of chemotaxis is complex, and existing models of chemotaxis are in general either mechanistic or phenomenological (Palsson and Othmer 2000; Hauser et al. 2007a; Figge et al. 2008). We use phenomenological models in this paper as our interest is in the feedback between receptors and an external chemokine field, and not so much in the detailed mechanism of chemotaxis. Phenomenological models are typically based on some variation of a persistent random walk biased in the direction of the chemokine gradient (Weiner 2002). To bridge between deterministic and stochastic motility models in continuous time, we use the classical Langevin process stochastic differential equation formalism for persistent random walks to model chemotaxis and reduce to a deterministic version by removing the Wiener noise component where appropriate. We also explicitly incorporate GPCR desensitization by an external chemokine field in the Langevin process model.

Chemokine receptors are regulated on multiple levels, and receptor dynamics can be complex (Lauffenburger and Linderman 1996). As an illustrative example (Beyer and Meyer-Hermann 2008) present a detailed formalism (comprising 6 differential equations with 13 free parameters) to model the dynamics of a single receptor type interacting with its chemokine. Implementing a model with that degree of complexity would focus attention on detail and detract from our intention to show that very simple mechanisms suffice to induce complex migratory behavior. We therefore chose to derive a new phenomenological model for the receptor incorporating just ligand binding, constitutive and binding-induced down-regulation, and de novo synthesis. The use of singular perturbation analysis leads to the formulation of a single equation to model the dynamics of each chemokine receptor.

Based on the considerations above, we propose a simple ODE model of individual B cells coupled to static chemokine fields. We used the model to investigate the range of dynamical behaviors exhibited in the presence of static chemokine field distributions representing the DZ and LZ of the germinal center. Our hypothesis was that study of a phenomenological model integrating spatial chemokine distributions, receptor regulation, and chemotaxis could provide a template for understanding the broad-stroke dynamics of B cells in the germinal center. This would complement the use of IBM simulations to fill in the fine details and reveal unexpected emergent behavior resulting from individual cell interactions. For our model, we chose to include just three components—a spatial distribution of chemokines in one dimension, a model for the regulation of chemokine receptors, and a chemotactic model for cell locomotion.

This manuscript describes the application of the spatially-driven ODE model to explore the migratory response of B cells to chemokines in the germinal center. We show that chemokine-induced receptor down-regulation and receptor-mediated chemotaxis in the presence of a simple fixed spatial distribution of the relevant chemokines is sufficient to induce complex migratory patterns, including intrazonal and interzonal oscillations.

Model Definition and Analysis
Static Chemokine Fields

The chemokine-driven ODE model is a deterministic nonlinear dynamical system in one spatial dimension, in which chemotaxis of a single cell is modulated by the levels of two chemokine receptors that are reciprocally regulated by the static 1D spatial distribution of their cognate chemokines.

Spatial Distribution of Chemokines

Chemokine fields are set up by the expression of chemokines by stromal cells in the germinal center with dynamics determined by diffusion, absorption and degradation, but in the steady state over short periods of time, we make the assumption that the chemokine field is static. We further simplify by assuming that each chemokine field has a Gaussian distribution, and only consider the dynamics along the axis that runs through both the follicle or germinal center centroids. Chemokines are modeled as functions of the cell displacement x—even though the field is static, cells with different displacements respond to the local chemokine field. This representation of chemokine fields as a function of cell displacement is flexible—it is possible to set up arbitrarily complex chemokine fields in this system if necessary to model in vivo measurements, for example, by using mixtures of Gaussians to represent multimodal fields.

Germinal Center Model

The chemokine concentration is given as a function of the cell displacement x. For the examples in the paper, we use Gaussian distributions f1 and f2 to represent the CXCL12 and CXCL13 chemokine distributions respectively, i.e.,

[Formula ID: Equ1]
[Formula ID: Equ2]
where ci and wi determine the height and width of the distribution, and k=(k1+k2)/2 is the half-distance between the dark and light zones in μm. The chemokine concentrations and gradients for CXCL12 and CXCL13 are shown in the first two panels in the top row of Fig. 1.

Toy Model for Receptor Regulation and Chemotaxis

We begin with a toy model for receptor regulation and chemotaxis to illustrate the basic requirements for chemokine-driven oscillations. We assume that the chemokine receptors are synthesized at a rate π and degraded at a rate δ. To couple the receptor dynamics to the chemokine field, we assume that receptors are also down-regulated at a rate proportional to the product of the receptor and the local cognate chemokine concentration. In other words, the chemokine drives the down-regulation of its receptor.

[Formula ID: Equ3]
[\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{dr_{i=1,2}}{dt} = \pi_{i} - r_{i} f_{i}(x) - \delta_{i} r_{i}$$\end{document}]

For chemotaxis, we assume that the cell velocity depends on the product of the receptor and the local cognate chemokine gradient, with a drag coefficient γ.

[Formula ID: Equ4]
[Formula ID: Equ5]

Analysis of Toy Model

In this section, we first analyze the toy model components individually to gain insight into the origin of specific dynamical behaviors, then integrate the components and explore the resulting system dynamics.

Regulation of Receptor Density

From Eq. (3), it follows that the steady state receptor density rSS at any given position x is given by

[Formula ID: Equ6]
[\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_i^{\text{SS}} = \frac{\pi_i}{\delta_i+f_i(x)}$$\end{document}]

In the rightmost upper panel of Fig. 1, we plot the steady state solution for the receptor density at a particular position. It is clear that the effect of binding-induced down-regulation is to decrease the receptor density the greatest where the chemokine concentration is highest. Where the level of cognate chemokine is low, synthesis of new receptor outpaces down-regulation, and the saturating density of receptor is achieved.

Stability of Cell Velocity

Next, we examine the dynamics of the velocity v with respect to relative changes in receptor concentration using a reduced undamped model (γ=0)

[Formula ID: Equ7]
[\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{dv}{dt} = s f_1'(x) + f_2'(x)$$\end{document}]
where we rescale so that s=r1/r2. We also assume that the chemokine fields f1(x) and f2(x) are standard normal distributions with centers set 1.5 units from the origin.

In the lower panels of Fig. 1, we plot the rate of change of the velocity as the position of a B cell is varied. There are three different sets of steady state solutions for v possible. When s, the ratio of the densities of the two chemokine receptors, is small (so that r2 dominates), there is a single steady state at the mean of f2(x). As s increases, a new steady state is created at the mean of f1(x) by a saddle-node bifurcation, and the system is bistable. As s continues to rise, the steady state at f2(x) vanishes in a reverse saddle node bifurcation, and the system becomes monostable again. This implies that under these conditions, the DZ, LZ or both can be equilibria for a B cell, depending only on the ratio of CXCR4 and CXCR5 expressed.

Coupling Receptor and Velocity Dynamics Results in Spontaneous Oscillations

Referring to the bottom panels of Fig. 1, we see how oscillations can arise from coupling of the receptor and velocity dynamics in the presence of opposing chemokine fields. Suppose a cell starts with a low density of CXCR4 and high CXCR5 at the CXCL12 peak. Under appropriate conditions, the only stable equilibrium is at the CXCL13 peak and the cell moves to the right (bottom left). When it reaches the CXCL13 peak, the chemokine drives the down-regulation of the CXCR5 receptor and CXCR4 is up-regulated. Now we are in the situation illustrated by the bottom right panel; the stable equilibrium at the CXCL13 peak vanishes, and the cell is forced to return to the CXCL12 peak, setting up a system where oscillations result.

Bifurcation Analysis

To better understand the conditions where the system exhibits oscillatory behavior, we can systematically study the dynamics under changes of parameters using software for continuation of equilibria (Dhooge et al. 2003; Clewley et al. 2007). Continuation software “follows” the equilibrium solution as some parameter is changed, and also checks for the occurrence of specific bifurcations at each parameter value. This allows us to identify parameter regions where interesting or desired system behavior is found, and provides insight into the parameter values where there is a qualitative change of behavior. Figure 2 shows the bifurcation plots from different parameter regions for the toy model.

Biologically Motivated Phenomenological Model for Receptor Regulation and Chemotaxis

The toy model described above shows that a combination of receptor adaptation (modeled as down-regulation) and receptor-mediated chemotaxis can give rise to autonomous oscillations in the presence of a suitable static chemokine field. By design, the model abstracts away all other biological considerations. In this section, we describe simple biologically-motivated models that accommodate standard mass-action kinetics for receptor dynamics and chemotaxis with saturable chemokine receptor signals, and show that these models preserve similar autonomous oscillations.

Model for GPCR Regulation

As discussed in a recent review by Bennett et al. (2011), the regulation of chemokine receptor levels is highly complex with multiple different processes that can affect GPCR levels and activity. However, the mechanism of migration is thought to be independent of transcription and regulated primarily by receptor trafficking dynamics in response to agonist binding (Schaeuble et al. 2012). Agonist-dependent desensitization in response to agonist binding results in GPCR endocytosis and degradation as discussed in the Introduction. Some of the internalized receptor may be recycled rather than degraded, and the path taken depends on both cell type and the duration of ligand engagement. While probably too slow a process to directly influence lymphocyte migration, new receptor synthesis is also essential for long-term maintenance of surface chemokine receptor levels. These processes of new receptor synthesis and agonist-induced internalization, recycling and degradation that together determine the dynamics of chemokine receptor expression are illustrated in Fig. 3. The mass-action kinetics corresponding to Fig. 3 are given by

[Formula ID: Equ8]
[Formula ID: Equ9]
[Formula ID: Equ10]
10
where the first-order degradation term τ for U is necessary to ensure that the unbound receptor remains finite in the absence of ligand. To simplify the model, we neglect the contribution of receptor recycling on the available intracellular pool. That is, we assume that πμB, and hence that I is constant. With these assumptions, we can derive the following model for the dynamics of the CXCR4 and CXCR5 receptors in the presence of their cognate chemokine ligands (full derivation given in Appendix A)
[Formula ID: Equ11]
11  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{dr_{i=1,2}}{dt} = \pi_{i} - \frac{r_{i} \tau_{i} \kappa_{i} f_{i}(x)}{1 + \kappa_{i} f_{i}(x)} - \delta_{i} r_{i}$$\end{document}]
where κ is a rescaled equilibrium association constant for GPCR:lignad interactions, and δ+τ is the removal rate for bound GPCR that incorporates the first-order degradation of unbound receptor U.

Chemotactic Model for Cell Locomotion

The model for B cell chemotaxis assumes that the cell velocity is governed by a chemotactic process biased by saturable chemokine receptor signals generated by receptor ligand interactions that depend on both ligand concentration and density. The velocity has a drag coefficient γ, and a tuning factor for the degree of responsiveness to the underlying chemokine field given by χ. The model equations (derived in Appendix B) are

[Formula ID: Equ12]
12
[Formula ID: Equ13]
13

We have set the effective equilibrium association constant ϵ in the chemotactic model to be distinct from the value κ in the receptor regulation model to allow for differential coupling of bound receptor to signal transduction pathways involved in the two processes.

Bifurcation Analysis of Phenomenological Model

The full model is reproduced below for convenience.

[Formula ID: Equ14]
14
[Formula ID: Equ15]
15
[Formula ID: Equ16]
16
[Formula ID: Equ17]
17

As with the toy model, the bifurcation analysis suggests that spontaneous oscillations only occur for a rather restricted set of parameter combinations. For example, the leftmost panel of Fig. 4 shows that spontaneous oscillations only occur when the distance separating the simulated DZ and LZ fall within a narrow range.

Diversity of Dynamical Behaviors in 1D

A surprisingly rich variety of periodic behavior can be found in the 1D ODE model system. The range of behaviors include direct passage to a steady state equilibrium, damped oscillations to steady state and a variety of oscillatory behaviors with periods ranging from minutes to many hours. These oscillations may be from DZ to LZ, within a single zone or have components of both small intrazonal and large inter-zonal circulations. Oscillations can be highly asymmetrical, with a disproportionate amount of time spent in a single zone. Oscillations may even be apparently chaotic. Representative examples of the numerical simulations of displacement over time illustrating the diversity of oscillations are shown in Fig. 5.

Individual-Based Model Simulations of Receptor Dynamics and Chemotaxis

Finally, we implemented the phenomenological model in a 3D IBM simulation of immune cells (Kepler and Chan 2007; Mitha et al. 2008), extending the chemotactic model to incorporate stochastic deviations. We show that very similar dynamical behavior is observed in the 3D simulation as in the simpler ODE models.

Stochastic Model for Chemotaxis

For the IBM, we rewrite the phenomenological model for chemotaxis as a stochastic differential equation in Ito form, giving rise to the following Langevin process equation

[Formula ID: Equ18]
18
[Formula ID: Equ19]
19
where dW is the differential Wiener process.

There are three main differences between the phenomenological model and IBM simulation—the IBM is in 3D while the phenomenological model is 1D; cells in the IBM have a stochastic chemotactic motility model rather than a deterministic one (i.e., σ≠0); and there are extrinsic forces in the IBM when cells collide with each other or environmental boundaries. While we have closed form solutions for the spatial distribution of chemokines in the two models shown here, the simulation system uses a spatially discretized numerical approximation in order to generalize to arbitrary (and potentially evolving) chemokine distributions. Numerically, the differences between the IBM simulation and 1D phenomenological models are the use of a three-dimensional grid to store chemokine concentrations and gradients (5 μm per side voxels with trilinear interpolation between voxel centroids) as compared with values given by the closed forms f1 and f2 in the ODE model. In addition, cells in the IBM can have more complex behaviors such as division, death, and activation and the possibility of collision-induced forces when multiple cells are simulated. In the IBM simulation, cells are also constrained to be within a specified volume.

Dynamical Behavior in Individual-Based Simulation

Figure 6 shows three snapshots of the 3D IBM simulation. Oscillatory behavior is preserved in the presence of external forces and stochasticity, although of course, the periodicity is no longer synchronized between cells.

In the absence of stochasticity (Fig. 7 middle panel), the 3D simulation model behavior corresponds very closely to that of the 1D phenomenological model (Fig. 7 top panel). With stochasticity (Fig. 7 bottom panel), the 3D simulation behavior begins to diverge. However, the dynamical analysis of the minimal model remains informative for the 3D simulation behavior.

Discussion

We have described a simple mathematical model of chemotaxis-driven B cell migration in the germinal center. The model incorporates a static chemokine field, chemokine-induced receptor modulation, and chemotaxis driven by the interaction of the chemokine receptor with the local chemokine concentration and gradient. The model is specified using coupled first-order differential equations, lending itself to detailed analysis using techniques from nonlinear dynamics. Using this basic setup, we investigated the dynamics of B cell migration under a simple chemokine field comprising of two Gaussian distributions representing CXCL12 in the light zone and CXCL13 in the dark zone of the germinal center.

In this simple germinal center model, we show that spontaneous oscillations between the light and dark zone can arise, and the periodicity can be tuned so that the residence times in the dark and light zones is consistent with experimental observations. An interesting prediction of the model is that for a fixed width of the chemokine fields, oscillations only occur for a narrow range of separations between the dark and light zone. When the light and dark zones are too close or far apart, no oscillations are observed. Oscillations can also be elicited in an alternative model where one receptor is fixed, and only one receptor is regulated by the chemokine field (supplementary Fig. A.1), but then the allowed range of separations is even narrower. This suggests that reciprocal regulation of both CXCR4 and CXCR5 receptors gives more robust oscillatory behavior than regulation of a single receptor.

While the simple mechanism of chemokine-driven receptor down-regulation is sufficient for inducing autonomous oscillations of some complexity, the extent to which such a mechanism contributes to B cell cycling in the germinal center is unknown. In fact, there is substantial evidence that B cell cycling in the germinal center is largely driven by extrinsic influences (e.g., B cell:FDC or B cell:T follicular helper cell interactions) that trigger differentiation programs regulating the expression of chemokine receptors. However, our model shows that surprisingly complex migratory patterns can emerge from very simple mechanisms, a recurring theme in the study of nonlinear dynamical systems. We believe that this provides a useful alternative perspective on the causal mechanisms of complex immune cell migration patterns, such as those observed in the germinal center.

This work was originally motivated by the desire to simplify IBM simulations of B cell behavior in order to gain insight into observed motility patterns and to facilitate parameter calibration. The 1D phenomenological model described in Eqs. (14)–(17) differ from the single cells in the 3D IBM simulations in the restriction to one dimension, the absence of a stochastic component, and the absence of collisions with other cells and the environment boundaries. However, we show that the phenomenological model effectively predicts the large-scale behavior of the IBM simulation when parameters are matched. Dynamical behaviors of interest can be rapidly identified in the phenomenological model configuration using bifurcation analysis and numerical simulations, and then studied in the more realistic 3D stochastic context with the IBM simulation using the same parameter values as the 1D phenomenological model. This is much more efficient than the brute-force search over parameter space otherwise necessary for IBM simulations, since such models are analytically intractable and highly demanding of computational resources. A caveat is that the extent to which such ODE-based model simplifications can replicate the dynamics of richer IBM that incorporate phenomena such as cell-cell interactions is not known. We conjecture that ODE models with mean-field approximations of cell-cell interactions will still be useful for providing insight into the parameters of these more challenging simulations and their calibration, and plan to investigate such approximations.

In conclusion, the ODE models for B cell motility described offer potential for a thorough analysis of the surprising complexity engendered by simple environment/cell interactions, and highlight the importance of considering the chemokine environment in understanding migration patterns of B cells. In addition, the ODE models provide flexibility to perform rapid prototyping of B cell migration dynamics, and may serve as a tractable bridge to more detailed IBM simulations.

Appendix A: Model for GPCR Regulation

Let r be the density of unbound GPCR, and r the density of GPCR bound by its ligand, whose concentration is l. The dynamical equations for this system are

[Formula ID: Equ20]
20
[Formula ID: Equ21]
21
where π is the rate of production of GPCR where, δ is the removal rate for GPCR in the unbound state, and δ+τ is the removal rate for the bound GPCR.
[Formula ID: Equa]
[\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\kappa}{\varepsilon}$$\end{document}]
is the forward rate constant for ligand binding,
[Formula ID: Equb]
[\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{\varepsilon}$$\end{document}]
is the rate constant for dissociation. We use these expressions to facilitate taking the limit where the binding and dissociation reactions are much faster than the cellular processes. We reduce the complexity of the dynamical system by taking some of its subprocesses as occurring on a much faster time-scale than others. These fast subprocesses are treated as if in equilibrium with the more slowly varying components, thus eliminating dynamic degrees of freedom. The mathematics used is singular perturbation theory (Jones 1995). In the present case, the justification for making the approximation comes from the measured rates for the subprocesses. The binding and dissociation of CXCL12 and CXCR4 has been characterized using surface plasmon resonance, giving kon=4.20×105 M s1, koff=8.24×10−3 s1, and KD=3.47×10−8 M, showing that the reactions equilibrate with a characteristic time of about 44 seconds for CXCL12 concentration equal to the reactions’ Kd (35 nM) (Vega et al. 2011). In contrast, the characteristic time for CXCR4 internalization upon binding by CXCL12 was found to be between 450 and 600 seconds (estimated from Signoret et al. 1997, Fig. 8A).

Define the total GPCR density rTr+r , whose rate equation is

[Formula ID: Equ22]
22  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{d{r_T}}}{dt} = \pi - \delta{r_T} - \tau{r^*}$$\end{document}]
Equation (21) can now be written
[Formula ID: Equ23]
23  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{d{r^*}}}{dt} = - ( {\delta + \tau} ){r^*} + \frac {\kappa}{\varepsilon} \bigl( {{r_T} - {r^*}} \bigr)l - \frac {1}{\varepsilon}{r^*}$$\end{document}]

We now perform the singular perturbation analysis, taking the limit ε→0. We first construct the “outer solution” by expanding the state variables as

[Formula ID: Equ24]
24  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}\begin{aligned}[c] {r_T} ( {t,\varepsilon} ) &\equiv{r_{T,0}} + \varepsilon {r_{T,1}} + O \bigl( {{\varepsilon^2}} \bigr) \hfill \\ {r^*} ( {t,\varepsilon} ) &\equiv r_0^* + \varepsilon r_1^* + O \bigl( {{\varepsilon^2}} \bigr) \end{aligned}\end{document}]
and then matching coefficients of ε in Eqs. (22) and (24). To lowest order, we have
[Formula ID: Equ25]
25  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}\begin{aligned}[c] \frac{{d{r_{T,0}}}}{dt} &= \pi - \delta{r_{T,0}} - \tau r_0^* \\ 0 &= \kappa \bigl( {{r_{T,0}} - r_0^*} \bigr)l - r_0^* \end{aligned}\end{document}]
the second of Eqs. (25) has solution
[Formula ID: Equ26]
26  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_0^* = \frac{{\kappa{r_{T,0}}l}}{{1 + \kappa l}}$$\end{document}]
so that the first of Eqs. (25) becomes
[Formula ID: Equ27]
27  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{d{r_{T,0}}}}{dt} = \pi - \biggl( {\delta + \frac{{\tau\kappa l}}{{1 + \kappa l}}} \biggr){r_{T,0}}$$\end{document}]

To impose initial conditions, we must compute the “inner solution” obtained by rescaling time as

[Formula ID: Equ28]
28  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T = \frac{t}{\varepsilon}$$\end{document}]
and letting
[Formula ID: Equ29]
29  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R_T} ( {T,\varepsilon} ) = {r_T} \bigl( {t(T, \varepsilon ),\varepsilon} \bigr)$$\end{document}]
etc.

Now, matching coefficients of ε in the resulting differential equations gives

[Formula ID: Equ30]
30  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}\begin{aligned}[c] \frac{{d{R_{T,0}}}}{dt} &= 0 \\ \frac{{dR_0^*}}{dt} &= \kappa \bigl( {{R_{T,0}} - R_0^*} \bigr)l - R_0^* \end{aligned}\end{document}]

Now the initial value problem in the inner solution has solution

[Formula ID: Equ31]
31  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R_{T,0}} ( T ) = {r_T} ( 0 )$$\end{document}]
and
[Formula ID: Equ32]
32  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0^* ( T ) = \biggl( {{r^*} ( 0 ) - \frac{{\kappa l{r_T}(0)}}{{1 + \kappa l}}} \biggr){e^{ - ( {1 + \kappa l} )T}} + \frac{{\kappa l{r_T}(0)}}{{1 + \kappa l}}$$\end{document}]

Matching the inner and outer solutions requires setting

[Formula ID: Equ33]
33  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r_{T,0}} ( 0 ) = \lim _{T \to\infty} {R_{T,0}} ( T ) = {r_{T,0}} ( 0 )$$\end{document}]
and
[Formula ID: Equ34]
34  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_0^* ( 0 ) = \lim_{T \to\infty} R_0^* ( T ) = \frac{{\kappa l{r_T}(0)}}{{1 + \kappa l}}$$\end{document}]

Equation (11) from the text is recovered by dropping subscripts and substituting the chemokine fields for the ligand concentration,

[Formula ID: Equ35]
35

Initial parameter values for the spatial measurements and receptor dynamics used for bifurcation analysis and numerical simulations were derived from the literature (Lin and Butcher 2008; Hauser et al. 2007a; Allen et al. 2004; Victora et al. 2010; Zigmond 1981; Hoffman et al. 1996; Ricart et al. 2011), or estimated when no experimental data was available.

Appendix B: Model for Cellular Locomotion

The model for cellular locomotion starts with a third-order Langevin process in Ito form:

[Formula ID: Equ36]
36
where X(t)∈ℝ3 is the cell’s position at time; V(t)∈ℝ3 and P(t)∈ℝ3 are the velocity and polarization, respectively. We use upper-case letters to remind us that these variables are stochastic. The effective drag coefficient is γ, and the polarization decorrelation rate is ζ. Φ is the external force exerted on the cell, and Γ is the signal due to an external orientation field. dW is a Wiener process generating fluctuations in the polarization, and σ controls the size of those fluctuations.

We proceed from here by providing a model for Γ in the case where the orientation field is due to an inhomogeneous chemokine distribution. We suppose that the binding of chemokine receptors on the cell’s surface generates a local signal, and the global orientation signal is the vector average of these signals over the cell’s surface.

[Formula ID: Equ37]
37  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varGamma\propto\int_{S}{{d^{2}}y \,{r^{*}} ( y )\widehat {n} ( y )}$$\end{document}]
where, as in Appendix A, r is the density of bound receptor, now considered a function of position y on the cell surface S, [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\widehat {n} ( y )$\end{document}] is the unit vector normal to the cell surface at surface point y. We use the same singular perturbation method to get
[Formula ID: Equ38]
38  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r^{*}} ( y )\propto\frac{f ( x+y )r}{1+\kappa f ( x+y )}=\frac{f ( x )r}{1+\kappa f ( x )}+r{y^{T}} \nabla\frac{f ( x )}{1+\kappa f ( x )}+O \bigl( {{\Vert y \Vert }^{2}} \bigr)$$\end{document}]
where x is the position of the center of the cell, κ is the equilibrium association constant, and f(x) is the concentration of chemokine at x. The expression to the right of the equals sign results from a Taylor expansion.

Substituting Eq. (38) into Eq. (37) gives

[Formula ID: Equ39]
39  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varGamma=\chi\frac{r\nabla f ( x )}{{{ [ 1+\kappa f ( x ) ]}^{2}}}$$\end{document}]
where the constant χ is the chemotactic coefficient.

Finally, we let the coefficients ζ,χ, and σ become very large while their ratios remain constant. We further assume that there are no external forces. In this limit, we get the system of equations

[Formula ID: Equ40]
40  [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}\begin{aligned}[c] X &=V\,dt \\ dV &=-\gamma V\,dt+\chi\frac{r\nabla f ( x )}{{{ [ 1+\kappa f ( x ) ]}^{2}}}\,dt+\sqrt{\gamma}\sigma \,dW \end{aligned}\end{document}]
where the parameters have been rescaled to give the form displayed.

Values for the motility parameters used in the simulation were calculated by fitting to data from 3D trajectory data of individual lymphocytes from 2-photon data (Kepler and Chan 2007).

Acknowledgements

This work was supported by NIH/NIAID research contract HHSN272201000053C (TB Kepler, PI).

Open Access