A Model for Migratory B Cell Oscillations from Receptor DownRegulation Induced by External Chemokine Fields.  
Jump to Full Text  
MedLine Citation:

PMID: 23296998 Owner: NLM Status: Publisher 
Abstract/OtherAbstract:

A longstanding 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 reentry 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, twophoton 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 postactivation germinal center. There is strong evidence that chemokines secreted by stromal cells and the regulation of cognate Gprotein 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 individualbased 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 individualbased 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 downregulation 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 cellimprinted poly(dimethylsiloxane) surface on the cellular activities o... 22992018  Metabolism of inorganic arsenic in intestinal epithelial cell lines. 23722198  Induction of apoptosis by highdose gold nanoparticles in nasopharyngeal carcinoma cells. 17065808  Activation of kinases upon volume changes: role in cellular homeostasis. 3383578  Quantitative effects of short and longterm vasectomy on mouse spermatogenesis and spe... 
Publication Detail:

Type: JOURNAL ARTICLE Date: 201318 
Journal Detail:

Title: Bulletin of mathematical biology Volume:  ISSN: 15229602 ISO Abbreviation: Bull. Math. Biol. Publication Date: 2013 Jan 
Date Detail:

Created Date: 201318 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:

Full Text  
Journal Information Journal ID (nlmta): Bull Math Biol Journal ID (isoabbrev): Bull. Math. Biol ISSN: 00928240 ISSN: 15229602 Publisher: SpringerVerlag, 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 pmcrelease 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/s1153801297999 
A Model for Migratory B Cell Oscillations from Receptor DownRegulation 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 
The evolution of highaffinity specific antibodies by longlived 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 longheld hypothesis of cyclic reentry 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}; MeyerHermann et al. ^{2001}). FDCs present antigen bound on Fc receptorcaptured 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 longlived plasma cells.
How B cells migrate in the lymph node is hence critical for understanding the generation of high affinity longlived 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 receptorligand interactions also modulate the B cell spatial distribution in the follicle (Gatto et al. ^{2011}). In the postactivation 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 Gprotein 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 GPCRmediated chemotaxis (Zigmond ^{1981}; Zigmond et al. ^{1982}). Estimated receptor levels are in the 10^{4} range (10,000–50,000) of receptors per cell (Zigmond ^{1981}). Upon ligand binding, GPCRs signal to Gprotein, become phosphorylated by GPCR kinase (effectively desensitizing the receptor by dissociating Gprotein subunits), and internalize via one of two major pathways. One internalization pathway is fast and involves clathrinmediated endocytosis. The other is slower and uses a lipidraft/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 twophoton 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, twophoton microscopy is restricted to the visualization of relatively small regions and short timespans. 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}; MeyerHermann et al. ^{2009}; Figge and MeyerHermann ^{2011}). As traditional ordinary differential equation (ODE) models ignore the spatial inhomogeneity of the chemokine fields, computational simulations of individualbased 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 MeyerHermann ^{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 bindinginduced downregulation, 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 broadstroke 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 spatiallydriven ODE model to explore the migratory response of B cells to chemokines in the germinal center. We show that chemokineinduced receptor downregulation and receptormediated 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.
The chemokinedriven 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.
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.
The chemokine concentration is given as a function of the cell displacement x. For the examples in the paper, we use Gaussian distributions f_{1} and f_{2} to represent the CXCL12 and CXCL13 chemokine distributions respectively, i.e.,
where c_{i} and w_{i} determine the height and width of the distribution, and k=(k_{1}+k_{2})/2 is the halfdistance 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.We begin with a toy model for receptor regulation and chemotaxis to illustrate the basic requirements for chemokinedriven 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 downregulated at a rate proportional to the product of the receptor and the local cognate chemokine concentration. In other words, the chemokine drives the downregulation of its receptor.
3 [
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 γ.
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.
From Eq. (3), it follows that the steady state receptor density r^{SS} at any given position x is given by
6 [
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 bindinginduced downregulation 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 downregulation, and the saturating density of receptor is achieved.
Next, we examine the dynamics of the velocity v with respect to relative changes in receptor concentration using a reduced undamped model (γ=0)
7 [
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 r_{2} dominates), there is a single steady state at the mean of f_{2}(x). As s increases, a new steady state is created at the mean of f_{1}(x) by a saddlenode bifurcation, and the system is bistable. As s continues to rise, the steady state at f_{2}(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.
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 downregulation of the CXCR5 receptor and CXCR4 is upregulated. 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.
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.
The toy model described above shows that a combination of receptor adaptation (modeled as downregulation) and receptormediated 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 biologicallymotivated models that accommodate standard massaction kinetics for receptor dynamics and chemotaxis with saturable chemokine receptor signals, and show that these models preserve similar autonomous oscillations.
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}). Agonistdependent 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 longterm maintenance of surface chemokine receptor levels. These processes of new receptor synthesis and agonistinduced internalization, recycling and degradation that together determine the dynamics of chemokine receptor expression are illustrated in Fig. 3. The massaction kinetics corresponding to Fig. 3 are given by
where the firstorder 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)11 [
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
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.
The full model is reproduced below for convenience.
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.
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 interzonal 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.
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.
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
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 threedimensional 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 f_{1} and f_{2} 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 collisioninduced forces when multiple cells are simulated. In the IBM simulation, cells are also constrained to be within a specified volume.
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.
We have described a simple mathematical model of chemotaxisdriven B cell migration in the germinal center. The model incorporates a static chemokine field, chemokineinduced 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 firstorder 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 chemokinedriven receptor downregulation 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 largescale 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 bruteforce 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 ODEbased model simplifications can replicate the dynamics of richer IBM that incorporate phenomena such as cellcell interactions is not known. We conjecture that ODE models with meanfield approximations of cellcell 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.
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
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.[
[
Define the total GPCR density r_{T}≡r+r^{∗} , whose rate equation is
22 [
23 [
We now perform the singular perturbation analysis, taking the limit ε→0. We first construct the “outer solution” by expanding the state variables as
24 [
25 [
26 [
27 [
To impose initial conditions, we must compute the “inner solution” obtained by rescaling time as
28 [
29 [
Now, matching coefficients of ε in the resulting differential equations gives
30 [
Now the initial value problem in the inner solution has solution
31 [
32 [
Matching the inner and outer solutions requires setting
33 [
34 [
Equation (11) from the text is recovered by dropping subscripts and substituting the chemokine fields for the ligand concentration,
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.
The model for cellular locomotion starts with a thirdorder Langevin process in Ito form:
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 uppercase 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.
37 [
38 [
Substituting Eq. (38) into Eq. (37) gives
39 [
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
40 [
Values for the motility parameters used in the simulation were calculated by fitting to data from 3D trajectory data of individual lymphocytes from 2photon data (Kepler and Chan ^{2007}).
This work was supported by NIH/NIAID research contract HHSN272201000053C (TB Kepler, PI).
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.. Germinalcenter 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.13652567.2011.03485.x21977995  
Beyer T.,MeyerHermann M.. Cell transmembrane receptors determine tissue pattern stabilityPhys. Rev. Lett.Year: 20081011410.1103/PhysRevLett.101.148102  
Bogle G.,Dunbar P. R.. Agentbased simulation of tcell 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., & MeyerHermann, M. (2011). Modelling intravital twophoton data of lymphocyte migration and interaction. Math. Models Immune Cell Biol., 121–139.  
Figge M. T.,Garin A.,Gunzer M.,KoscoVilbois M.,Toellner K. M.,MeyerHermann M.. Deriving a germinal center lymphocyte migration model from twophoton 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.,MeierSchellersheim M.,NitaLazar A.,Fraser I. D. C.. Systems biology in immunology—a computational modeling perspectiveAnnu. Rev. Immunol.Year: 20112952710.1146/annurevimmunol03040910131721219182  
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 Bcell helpBloodYear: 200510661924193110.1182/blood200411449415899919  
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 germinalcenter 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 germinalcentre developmentNat. Rev. Immunol.Year: 20077749950410.1038/nri212017589541  
Hoffman J. F.,Linderman J. J.,Omann G. M.. Receptor upregulation, 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 reentry of germinal center B cells and the efficiency of affinity maturationImmunol. TodayYear: 199314841241510.1016/01675699(93)90145B8397781  
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  
MeyerHermann M.,Deutsch A.,OrGuil M.. Recycling probability and dynamical properties of germinal center reactionsJ. Theor. Biol.Year: 2001210326528510.1006/jtbi.2001.229711397129  
MeyerHermann M.,Figge M. T.,Toellner K. M.. Germinal centres seen through the mathematical eye: Bcell 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 cellbased 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 Bcell 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.,KoscoVilbois 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.1600065X.2012.01124.x22500831  
Signoret N.,Oldridge J.,PelchenMatthews A.,Klasse P. J.,Tran T.,Brass L. F.,Rosenkilde M. M.,Schwartz T. W.,Holmes W.,Dallas W.,et al. Phorbol esters and SDF1 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íguezFrade J. M.,Calle A.,RodríguezFernández J. L.,Lechuga L. M.,Rodríguez J. F.,GutiérrezGallego R.,et al. Technical advance: surface plasmon resonancebased 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.,MeyerHermann 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/S09550674(02)00310111891119  
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=r_{1}/r_{2} 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 f_{1}(x) and f_{2}(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,r_{1}=1,r_{2}=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 volumerendered, 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:

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