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

This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.


References
Allen C. D. C.,Ansel K. M.,Low C.,Lesley R.,Tamamura H.,Fujii N.,Cyster J. G.. Germinal center dark and light zone organization is mediated by CXCR4 and CXCR5Nat. Immunol.Year: 20045994395210.1038/ni110015300245
Allen C. D. C.,Okada T.,Cyster J. G.. Germinal-center organization and cellular dynamicsImmunityYear: 200727219020210.1016/j.immuni.2007.07.00917723214
Allen C. D. C.,Okada T.,Tang H. L.,Cyster J. G.. Imaging of germinal center selection events during affinity maturationScienceYear: 2007315581152810.1126/science.113673617185562
Beltman J. B.,Allen C. D. C.,Cyster J. G.,De Boer R. J.. B cells within germinal centers migrate preferentially from dark to light zoneProc. Natl. Acad. Sci.Year: 201110821875510.1073/pnas.110155410821555569
Bennett L. D.,Fox J. M.,Signoret N.. Mechanisms regulating chemokine receptor activityImmunologyYear: 2011134324625610.1111/j.1365-2567.2011.03485.x21977995
Beyer T.,Meyer-Hermann M.. Cell transmembrane receptors determine tissue pattern stabilityPhys. Rev. Lett.Year: 20081011410.1103/PhysRevLett.101.148102
Bogle G.,Dunbar P. R.. Agent-based simulation of t-cell activation and proliferation within a lymph nodeImmunol. Cell Biol.Year: 200988217217910.1038/icb.2009.7819884904
Clewley, R. H., Sherwood, W. E., LaMar, M. D., & Guckenheimer, J. M. (2007). PyDSTool, a software environment for dynamical systems modeling. http://pydstool.sourceforge.net.
Dhooge A.,Govaerts W.,Kuznetsov Y. A.. Matcont: a MATLAB package for numerical bifurcation analysis of ODEsACM Trans. Math. Softw.Year: 200329214116420008801070.6557410.1145/779359.779362
Figge M. T.. Stochastic discrete event simulation of germinal center reactionsPhys. Rev. EYear: 2005715217182610.1103/PhysRevE.71.051907
Figge, M. T., & Meyer-Hermann, M. (2011). Modelling intravital two-photon data of lymphocyte migration and interaction. Math. Models Immune Cell Biol., 121–139.
Figge M. T.,Garin A.,Gunzer M.,Kosco-Vilbois M.,Toellner K. M.,Meyer-Hermann M.. Deriving a germinal center lymphocyte migration model from two-photon dataJ. Exp. Med.Year: 2008205133019302910.1084/jem.2008116019047437
Gatto D.,Wood K.,Brink R.. Ebi2 operates independently of but in cooperation with CXCR5 and CCR7 to direct b cell migration and organization in follicles and the germinal centerJ. Immunol.Year: 201118794621462810.4049/jimmunol.110154221948984
Germain R. N.,Meier-Schellersheim M.,Nita-Lazar A.,Fraser I. D. C.. Systems biology in immunology—a computational modeling perspectiveAnnu. Rev. Immunol.Year: 20112952710.1146/annurev-immunol-030409-10131721219182
Hardtke S.,Ohl L.,Förster R.. Balanced expression of CXCR5 and CCR7 on follicular T helper cells determines their transient positioning to lymph node follicles and is essential for efficient B-cell helpBloodYear: 200510661924193110.1182/blood-2004-11-449415899919
Hauser A. E.,Junt T.,Mempel T. R.,Sneddon M. W.,Kleinstein S. H.,Henrickson S. E.,von Andrian U. H.,Shlomchik M. J.,Haberman A. M.. Definition of germinal-center B cell migration in vivo reveals predominant intrazonal circulation patternsImmunityYear: 200726565566710.1016/j.immuni.2007.04.00817509908
Hauser A. E.,Shlomchik M. J.,Haberman A. M.. In vivo imaging studies shed light on germinal-centre developmentNat. Rev. Immunol.Year: 20077749950410.1038/nri212017589541
Hoffman J. F.,Linderman J. J.,Omann G. M.. Receptor up-regulation, internalization, and interconverting receptor statesJ. Biol. Chem.Year: 199627131183941840410.1074/jbc.271.31.183948702483
Jones, C. (1995). Geometric singular perturbation theory. Dyn. Syst., 44–118.
Kepler T. B.,Chan C.. Spatiotemporal programming of a simple inflammatory processImmunol. Rev.Year: 2007216115316317367341
Kepler T. B.,Perelson A. S.. Cyclic re-entry of germinal center B cells and the efficiency of affinity maturationImmunol. TodayYear: 199314841241510.1016/0167-5699(93)90145-B8397781
Kepler T. B.,Perelson A. S.,et al. Somatic hypermutation in B cells: an optimal control treatmentJ. Theor. Biol.Year: 19931641376410.1006/jtbi.1993.11398264243
Kleinstein S. H.. Toward quantitative models of germinal center dynamicsYear: 2002PrincetonPrinceton University
Lauffenburger D. A.,Linderman J. J.. Receptors: models for binding, trafficking, and signalingYear: 1996LondonOxford University Press
Lin F.,Butcher E. C.. Modeling the role of homologous receptor desensitization in cell gradient sensingJ. Immunol.Year: 2008181128335834319050250
Meyer-Hermann M.,Deutsch A.,Or-Guil M.. Recycling probability and dynamical properties of germinal center reactionsJ. Theor. Biol.Year: 2001210326528510.1006/jtbi.2001.229711397129
Meyer-Hermann M.,Figge M. T.,Toellner K. M.. Germinal centres seen through the mathematical eye: B-cell models on the catwalkTrends Immunol.Year: 200930415716410.1016/j.it.2009.01.00519282244
Mitha, F., Lucas, T., Feng, F., Kepler, T. B., & Chan, C. (2008). The multiscale systems immunology project: software for cell-based immunological simulation. Source Code Biol. Med., 3.
Palsson E.,Othmer H. G.. A model for individual and collective cell movement in dictyostelium discoideumProc. Natl. Acad. Sci.Year: 200097191044810.1073/pnas.97.19.1044810984537
Reif K.,Ekland E. H.,Ohl L.,Nakano H.,Lipp M.,Förster R.,Cyster J. G.. Balanced responsiveness to chemoattractants from adjacent zones determines B-cell positionNatureYear: 20024166876949910.1038/416094a11882900
Ricart B. G.,John B.,Lee D.,Hunter C. A.,Hammer D. A.. Dendritic cells distinguish individual chemokine signals through CCR7 and CXCR4J. Immunol.Year: 201118615310.4049/jimmunol.100235821106854
Schaeuble K.,Hauser M. A.,Rippl A. V.,Bruderer R.,Otero C.,Groettrup M.,Legler D. F.. Ubiquitylation of the chemokine receptor CCR7 enables efficient receptor recycling and cell migrationJ. Cell Sci.Year: 20121254463447410.1242/jcs.09751922797918
Schwickert T. A.,Lindquist R. L.,Shakhar G.,Livshits G.,Skokos D.,Kosco-Vilbois M. H.,Dustin M. L.,Nussenzweig M. C.. In vivo imaging of germinal centres reveals a dynamic open structureNatureYear: 20074467131838710.1038/nature0557317268470
Shlomchik M. J.,Weisel F.. Germinal center selection and the development of memory B and plasma cellsImmunol. Rev.Year: 20122471526310.1111/j.1600-065X.2012.01124.x22500831
Signoret N.,Oldridge J.,Pelchen-Matthews A.,Klasse P. J.,Tran T.,Brass L. F.,Rosenkilde M. M.,Schwartz T. W.,Holmes W.,Dallas W.,et al. Phorbol esters and SDF-1 induce rapid endocytosis and down modulation of the chemokine receptor CXCR4J. Cell Biol.Year: 1997139365166410.1083/jcb.139.3.6519348282
Vega B.,Muñoz L. M.,Holgado B. L.,Lucas P.,Rodríguez-Frade J. M.,Calle A.,Rodríguez-Fernández J. L.,Lechuga L. M.,Rodríguez J. F.,Gutiérrez-Gallego R.,et al. Technical advance: surface plasmon resonance-based analysis of CXCL12 binding using immobilized lentiviral particlesJ. Leukoc. Biol.Year: 201190239940810.1189/jlb.101056521593136
Victora G. D.,Schwickert T. A.,Fooksman D. R.,Kamphorst A. O.,Meyer-Hermann M.,Dustin M. L.,Nussenzweig M. C.. Germinal center dynamics revealed by multiphoton microscopy with a photoactivatable fluorescent reporterCellYear: 2010143459260510.1016/j.cell.2010.10.03221074050
Weiner O. D.. Regulation of cell polarity during eukaryotic chemotaxis: the chemotactic compassCurr. Opin. Cell Biol.Year: 2002142196202260092410.1016/S0955-0674(02)00310-111891119
Zigmond S. H.. Consequences of chemotactic peptide receptor modulation for leukocyte orientationJ. Cell Biol.Year: 198188364464710.1083/jcb.88.3.6446260816
Zigmond S. H.,Sullivan S. J.,Lauffenburger D. A.. Kinetic analysis of chemotactic peptide receptor modulationJ. Cell Biol.Year: 1982921344310.1083/jcb.92.1.346276415

Figures

[Figure ID: Fig1]
Fig. 1 

Characterization of static chemokine fields and dynamics of toy model. The upper left panel shows the static concentration of CXCL12 and CXCL13, and the upper center panel shows the gradient of CXC12 and CXC13. The upper right panel shows the steady state solutions for the receptor density induced by local chemokine concentration, with dashed curves representing chemokine concentrations (blue = CXCL12, red = CXCL13) and solid curves representing steady state receptor densities if a cell was kept fixed at that position (blue = CDCR4, red = CXCR5). The lower panels show the steady state solutions for the velocity v as s=r1/r2 varies. The black curve shows the rate of change of v, and the red circles indicate stable steady states. In the left panel, s=0 and the system only has a single steady state at x≈1.5. In the middle panel, s=1 and the system is bistable. Finally, in the right panel, s=4.5 and the system has lost the steady state at x≈1.5 and only the steady state at x≈−1.5 remains. Scaled chemokine concentrations represented by f1(x) and f2(x) are shown as dashed lines for reference (Color figure online)



[Figure ID: Fig2]
Fig. 2 

Bifurcation plots and sample trajectories for the toy model. Left panel shows the 1D bifurcation plot as the separation k between DZ and LZ is varied. Oscillations arise as k is increased above the threshold labeled as H1 where a supercritical Hopf bifurcation occurs. The dotted line segment indicates the parameter range for k where the equilibrium solution is unstable and oscillations exist. At approximately k=39.5, a subcritical Hopf bifurcation occurs at H2. The equilibrium becomes stable again at the saddle node bifurcation LP2 and there are no oscillations beyond this value. It is also possible to go beyond single parameters and identify parameter combinations where bifurcations occur as illustrated by the middle panel that shows the continuation of the Hopf bifurcation as k and the rate of synthesis of CXCR4 (π1) are simultaneously varied. The right panel shows sample trajectories for the pairs of (k, π1) parameters indicated by crosses in the middle panel, with trajectories for each parameter pair plotted matching the same color cross. Fixed parameters are chosen to be within reasonable biological ranges and to generate oscillations on the correct time scale (hours): c=20, w=15.7, π2=1, δ1=0.001, δ2=0.001, γ=1. In the left panel, π1=0.3, while in the center panel, the red (k=40, π1=0.35), green (k=45, π1=0.58) and blue (k=50, π1=0.805) crosses represent (k,π1) pairs that have oscillatory behaviors (Color figure online)



[Figure ID: Fig3]
Fig. 3 

Schematic showing model for regulation of chemokine receptor levels on the cell surface. Unbound receptors (U) and bound receptors (B) levels are determined by association or dissociation with chemokine (shown as red ellipse) with rates k+ and k, respectively. Bound receptors are endocytosed to become internalized receptors (I) with rate μ and internalized receptors can either be degraded with rate δ or recycled with rate β. Finally, new receptor synthesis occurs with rate π (Color figure online)



[Figure ID: Fig4]
Fig. 4 

Bifurcation plots and sample trajectories for the reduced germinal center model. Left panel shows the 1D bifurcation plot as the separation k between DZ and LZ is varied. Oscillations arise as k is increased above the threshold labeled as H1 where a supercritical Hopf bifurcation occurs. The dotted line segment indicates the parameter range for k where the equilibrium solution is unstable and oscillations exist. At approximately k=52, a subcritical Hopf bifurcation occurs at H2. The equilibrium becomes stable again at the saddle node bifurcation LP2 and there are no oscillations beyond this value. It is also possible to go beyond single parameters and identify parameter combinations where bifurcations occur as illustrated by the right panel that shows the continuation of the Hopf bifurcation as k and the rate of synthesis of CXCR4 (π1) are simultaneously varied. Fixed parameters: c=10, w=25, π2=0.1, τ1=0.06, τ2=0.06, δ1=0.006, δ2=0.006, κ1=1, κ2=0.1, ϵ1=0.3, ϵ2=0.3, χ=28, ζ=1, γ=5. In the left panel, π1=0.15, while in the right panel, the cyan (k=43, π1=0.9), black (k=44, π1=0.75), and green (k=50, π1=0.15) circles represent (k, π1) pairs that have oscillatory behavior, while the magenta (k=44, π1=0.9), blue (k=50, π1=0.25), and red (k=50, π1=0.05) circles represent (k, π1) pairs that have steady state solutions. In particular, the green circle (k=50, π1=0.15) corresponds to a solution in the unstable equilibrium region of the left panel, where oscillatory behavior is predicted. Sample trajectories corresponding to these parameter pairs are shown in Fig. 7 (Color figure online)



[Figure ID: Fig5]
Fig. 5 

Diversity of oscillations with the phenomenological model. In each case, initial conditions were x=0,v=1,r1=1,r2=1. Top row: Left panel shows symmetrical oscillations, center panel shows asymmetrical oscillations and right panel shows oscillations confined to one zone. Bottom row: Left panel shows “nested” large and small oscillations in both zones, center panel shows “nested” small oscillations only in one zone, and right panel shows chaotic oscillations. Parameters: Top left (c=10, w=25, π1=0.15, π2=0.15, τ1=0.06, τ2=0.06, δ1=0.006, δ2=0.006, κ1=0.5, κ2=0.5, ϵ1=0.3, ϵ2=0.3, χ=28, γ=5, k=50); top center (c=10, w=25, π1=0.23, π2=0.15, τ1=0.06, τ2=0.06, δ1=0.006, δ2=0.006, κ1=0.5, κ2=0.5, ϵ1=0.3, ϵ2=0.3, χ=28, γ=5, k=50); top right (c=10, w=25, π1=0.15, π2=0.15, τ1=0.06, τ2=0.06, δ1=0.006, δ2=0.006, κ1=3.6, κ2=0.5, ϵ1=0.3, ϵ2=0.3, χ=28, γ=5, k=50); bottom left (c=5, w=22, π1=0.02, π2=0.004, τ1=0.005, τ2=0.01, δ1=0.001, δ2=0.001, κ1=10, κ2=10, ϵ1=3, ϵ2=0.2, χ=26.67, γ=0.03, k=45); bottom center (c=5, w=22, π1=0.05, π2=0.006, τ1=0.04, τ2=0.015, δ1=0.001, δ2=0.001, κ1=10, κ2=10, ϵ1=6, ϵ2=0.1, χ=3.33, γ=0.03, k=45); bottom right (c=5, w=22, π1=0.05, π2=0.0056, τ1=0.005, τ2=0.018, δ1=0.001, δ2=0.001, κ1=5.5, κ2=10, ϵ1=6, ϵ2=0.1, χ=16.17, γ=0.03, k=45)



[Figure ID: Fig6]
Fig. 6 

Snapshots of IBM GC B cell simulation. All B cells were started in the lower zone (DZ) and undergo spontaneous periodic motion between DZ and LZ as predicted by the phenomenological model. Chemokine concentrations in the DZ and LZ are volume-rendered, and the cell color shows density of CXCR4 receptors as indicated by the color bar (Color figure online)



[Figure ID: Fig7]
Fig. 7 

Comparison of sample oscillations from the 1D phenomenological model (top), IBM simulation with no stochasticity (middle) and IBM simulation with stochasticity (bottom). Top—Sample trajectories for the six (k, π1) parameter pairs shown in Fig. 4 found using a numerical ODE solver. Middle—Sample trajectories for the six (k, π1) parameter pairs shown in 4 obtained by running the 3D simulation for a single cell using identical parameter values ″ with σ=0. There is excellent agreement with the ODE solutions shown above. Very minor differences (e.g., slight change in periodicity) are likely due to spatial discretization (5 μm per side voxels with trilinear interpolation) and lower temporal resolution of the 3D simulation compared with the adaptive numerical integrator used in the ODE solution. Bottom—Sample trajectories for the six (k, π1) parameter pairs shown in Fig. 4 obtained by running the 3D simulation for a single cell again with identical parameters except for σ=2. The overall qualitative behavior is similar to that seen with σ=0, but the presence of stochasticity reveals the nearby attractor structure (see sporadic spikes in top row, middle column, and bottom row, last column) reminiscent of stochastic resonance



Article Categories:
  • Original Article


Previous Document:  Coarse Grained Normal Mode Analysis vs. Refined Gaussian Network Model for Protein Residue-Level Str...
Next Document:  Calculation of the Outcomes of Remodeling of Arteries Subjected to Sustained Hypertension Using a 3D...