Dimensionality of Motion and Binding Valency Govern Receptor–Ligand Kinetics As Revealed by Agent-Based Modeling

Mathematical modeling and computer simulations have become an integral part of modern biological research. The strength of theoretical approaches is in the simplification of complex biological systems. We here consider the general problem of receptor–ligand binding in the context of antibody–antigen binding. On the one hand, we establish a quantitative mapping between macroscopic binding rates of a deterministic differential equation model and their microscopic equivalents as obtained from simulating the spatiotemporal binding kinetics by stochastic agent-based models. On the other hand, we investigate the impact of various properties of B cell-derived receptors—such as their dimensionality of motion, morphology, and binding valency—on the receptor–ligand binding kinetics. To this end, we implemented an algorithm that simulates antigen binding by B cell-derived receptors with a Y-shaped morphology that can move in different dimensionalities, i.e., either as membrane-anchored receptors or as soluble receptors. The mapping of the macroscopic and microscopic binding rates allowed us to quantitatively compare different agent-based model variants for the different types of B cell-derived receptors. Our results indicate that the dimensionality of motion governs the binding kinetics and that this predominant impact is quantitatively compensated by the bivalency of these receptors.

1. Introduction

In recent decades, computational biology has developed into an autonomous scientific discipline that has become indispensable for contemporary biological research. Major contributions of computational biology comprise: (i) directing studies by providing insights that cannot otherwise be obtained in wet-lab experiments, (ii) advancing biological research toward a quantitative science through large-scale computations, and (iii) generating experimentally testable hypotheses through simulations of mathematical models.

The strength of mathematical modeling is actually in the simplification of complex processes by focusing on the most relevant aspects of a system. The art of modeling is in the appropriate choice of a mathematical approach that describes all existing experimental data and still can make relevant predictions. At this point a reasonable compromise has to be made between the level of system complexity that is transferred into the mathematical model and the feasibility of simulations with regard to computational resources.

Models based on ordinary differential equations (ODE) are presumably most frequently applied in biological research, even though this modeling approach is only valid if the system under consideration consists of large amounts of constituents, e.g., molecules, that are homogeneously distributed or well stirred in some spatial environment (1). This is because ODE models do not explicitly account for any spatial aspects of a system and changes in system variables, e.g., concentrations of molecules, are consequently described by functions of time that are continuous and deterministic. However, these assumptions, which may be typically appropriate for chemical systems, are for biological systems at best applicable from a macroscopic point of view. In these macroscopic models the biological processes are characterized by two specific types of parameters, which are referred to as rates or reaction rates. Rates characterize unimolecular processes that occur spontaneously and have unit 1/time. Reactions involving two types of molecules, i.e., bimolecular processes, are characterized by reaction rates with unit 1/(concentration × time). Typical experimental assays to determine these macroscopic rates for uni- and bimolecular processes are the adhesion frequency assay and the surface plasmon resonance assay (2). The advantage of ODE models is that they are based on a minimal set of parameters and can be formulated with relative ease (1, 3), which makes them belonging to the so-called simple modeling approaches (4). Deterministic ODE models may be extended to account for the stochasticity of chemical reactions in solution. Various numerical schemes have been introduced by Gillespie to sample the underlying master equation for the probability to find the system in a particular state at a given time (5). These are referred to as the direct method (6) and the first reaction method (7) and were later advanced for computational speed-up with the next reaction method by Gibson and Bruck (8). Albeit more detailed than deterministic ODE models, all these approaches have in common that a macroscopic viewpoint on the system is taken.

In contrast, agent-based models (ABMs), which belong to the so-called detailed modeling approaches (4), consider biological systems from a microscopic viewpoint by taking details of their individual constituents in space and time into account. A system’s constituents, e.g., molecules and/or cells, are represented by agents in the model and their motion in a specific spatial environment as well as their stochastic interactions with other agents are monitored in the simulations. In this microscopic modeling approach, all reactions are performed with a specific probability per time-step. This implies that not only the rates for unimolecular processes are measured in unit 1/time, but also the reaction rates for bimolecular processes, because the microscopic reactions are between two single molecules and not between concentrations of molecules as is the case for macroscopic ODE models. The microscopic rates for molecular interactions could be experimentally measured using thermal fluctuation assays (2). However, the level of detail represented by ABM comes at the price of a relatively large number of model parameters, which may be unknown and/or even inaccessible to experiment (1, 9), and simulations of ABM are typically associated with a high computational load (10, 11).

In this study, we focus on specific receptor–ligand (RL) binding, i.e., antibody–antigen binding as a central part of the adaptive immune response, and model this process in a comparative fashion by ODE models and by ABM. Binding between receptors and ligands represents an essential process in the immune system by which important information is transferred. For example, in the process termed opsonization, pathogen-derived antigens can be neutralized and labeled by antibodies for removal from the organism. Antibodies are soluble molecules that play a key role in the humoral response of adaptive immunity (12), because they can bind antigens with high affinity and can provide life-long protection against specific antigens. Of interest, antibodies do also exist as membrane-anchored molecules on B lymphocytes and are then referred to as B cell receptors (BCR). Binding of cognate antigen by BCR activates naïve B cells in lymphoid organs, such as spleen and lymph node (13), and this may initiate a germinal center (GC) reaction for antibody affinity maturation (12). During a GC reaction, B cells are proliferating and mutating their BCR followed by the selection of B cells with BCR that have high affinities to presented antigens. B cells with BCR that successfully accomplished the selection procedure differentiate into plasma cells that produce large amounts of these BCR as soluble antibodies. The GC reaction has been the subject of various interdisciplinary studies combining experimental and theoretical investigations (5, 14–17). In particular, it could be shown that the GC reaction is not only initiated by antigen binding to BCR on B cells, but that its termination is as well regulated by the high-affinity antibodies produced in soluble form (18). Taken together, antibodies represent a prime example for this study because of three reasons: (i) they exist as soluble as well as membrane-anchored receptors, (ii) they have a peculiar Y-shaped morphology that raises the question on its impact on RL binding as compared to spherically shaped receptors, and (iii) they have two binding sites and can bind antigen mono- or bivalently. The computational biology approach that is pursued in this study allows investigating the relative importance of receptor morphology, binding valency and dimensionality of motion that depends on receptors being soluble or membrane anchored on a cell. Applying different modeling approaches, e.g., ODE models and ABM, in a comparative fashion enables a quantitative mapping of the macroscopic and microscopic viewpoint on RL binding dynamics.

2. Materials and Methods

2.1. Microscopic Modeling of Receptor–Ligand Binding

Agent-based models (ABMs) are widely used in computational biology to simulate processes at the microscopic scale (9–11, 19). The individual constituents of the biological system under consideration are represented as agents that can move in a defined spatial environment and can interact with each other according to specific rules. We studied receptor–ligand (RL) binding and, in particular, the impact of specific receptor properties on the dynamics of the binding process. While ligands were modeled as molecules in solution with spherical shape, we considered receptors with different morphologies, i.e., being either spherically shaped (O) or Y-shaped (Y), and in settings with different dimensionality of motion, i.e., receptors in solution (SOL) or membrane anchored (MEM) on the surface of a cell. The four combinations of receptor properties are depicted in Figures 1 and 2, and give rise to four different ABM variants. These are denoted by their receptor properties, respectively, as O-SOL (see Figures 1A and 2A), O-MEM (see Figures 1B and 2B), Y-SOL (see Figures 1C and 2C), and Y-MEM (see Figures 1D and 2D). Simulations of the different ABM variants are shown in Videos S1–S5 in Supplementary Material. While in what follows we describe the general setup of the ABM, a detailed overview of the model parameters and of their corresponding values is provided in the Table S1 in Supplementary Material.

FIGURE 1

Figure 1. Schemes of ABM variants for receptor–ligand binding. The ABM variants are composed of the same spherical environment (large gray sphere) containing a spherical cell (small gray sphere) at the center. Ligands (orange) are always soluble, whereas receptors (blue) are studied in the variants: spherical receptor morphology in (A) soluble (O-SOL) or (B) membrane-anchored (O-MEM) form and Y-shaped receptor morphology in (C) soluble (Y-SOL) or (D) membrane-anchored (Y-MEM) form.

FIGURE 2

Figure 2. Detailed representation of receptor morphologies. Two-dimensional projection of three-dimensional receptors in ABM variants (A) O-SOL, (B) O-MEM, (C) Y-SOL, and (D) Y-MEM. Each receptor consists of binding spheres and collision spheres that may be overlapping in position and size. Ligands can bind by encountering a receptor’s binding sphere but are prohibited to penetrate receptors by the collision spheres. Details on the parameter values are provided in Table S1 in Supplementary Material.

2.1.1. Model System

In this study, we considered the model system of a B cell with Y-shaped B cell receptors (BCR), because these receptors do as well exist in a soluble form as antibodies. In the ABM, BCR with their Fab-fragments as binding sites are represented by a cylindrical stem with two cylindrical arms and spherical binding regions at the distal sides, which are hereafter referred to as binding spheres. A schematic representation of the BCR in soluble and membrane-anchored form is shown, respectively, in Figure 2C for ABM variant Y-SOL and in Figure 2D for ABM variant Y-MEM. The binding spheres on top of each arm represent the active binding sites of the BCR, whose surface areas are estimated from the size of Fv-regions, i.e., the variable parts of the BCR Fab-arms. Thus, the binding spheres implicitly account for the attractive short-range interactions between the binding sites of receptors and ligands (20–23). For the reason of comparison between BCR and spherically shaped receptors, we set the values of binding radii such that the effective area of all binding spheres are of comparable size, as can be inferred from the relative receptor sizes in Figures 2A,B for ABM variants O-SOL and O-MEM, respectively. For the same reason, when comparing Y-shaped and spherically shaped receptors, we impose the condition that receptors can only bind one ligand at a time. In addition, we also compared Y-shaped receptors that can bind mono- and bivalently.

2.1.2. Molecular Diffusion and Interaction

Receptors and ligands perform diffusive motion in the ABM. The corresponding diffusion coefficients can vary by orders of magnitude for soluble and membrane-anchored receptors. Diffusion coefficients were estimated based on the Stokes-Einstein equation (24) and the values for the corresponding ABM variants (see Table S1 in Supplementary Material) were calculated as outlined in Supplementary Material. In this study, we aim to investigate the impact of the dimensionality of motion for different receptor morphologies during the process of RL binding. In the ABM, molecules with diffusion coefficient D move per time step Δt the specific distance $\mathrm{\Delta }s=\sqrt{2dD\mathrm{\Delta }t}$ in a direction of the d-dimensional space that is chosen from a uniformly random distribution. This motion involves also a random rotation of Y-shaped receptors around their two axes in a spherically uniform fashion.

Two types of interaction processes are possible in the ABM: binding of receptor and ligand to form a molecular complex and dissociation of such a complex into individual receptor and ligand. The latter process occurs with rate $k_{\mathit{off}}^{micro}$ and translates into the probability $p_{\mathit{off}}^{micro}=k_{\mathit{off}}^{micro}\mathrm{\Delta }t$ that a complex dissociates during one time step Δt. In this study, we set the microscopic and macroscopic dissociation rates to be equal, i.e., $k_{\mathit{off}}^{micro}=k_{\mathit{off}}^{macro}$. As analyzed in detail in Supplementary Material, this approach is valid for typical parameter values of antibody–antigen dissociation rates, implying that dissociation and rebinding are relatively rare processes. On the other hand, binding of diffusing receptor and ligand requires that these molecules first encounter each other in the spatial environment. Then, upon contact of the ligand with the respective binding sphere of a receptor, binding occurs with probability $p_{on}^{micro}=k_{on}^{micro}\mathrm{\Delta }t$, where $k_{on}^{micro}$ denotes the microscopic binding rate with unit s⁻¹. Note, that this rate is conceptually different from the macroscopic reaction rate $k_{on}^{macro}$ with unit µm³ s⁻¹, because the latter incorporates the process of encounter of molecules in a spatially homogeneous system by their concentrations. In this study, we establish a relation between $k_{on}^{micro}$ and $k_{on}^{macro}$ by mapping the microscopic and macroscopic RL binding kinetics onto each other.

2.1.3. Implementation and Simulation

We implemented the ABM in a spherical environment with the cell positioned at its center and for reasons of comparison this was the same in all four ABM variants. The boundary condition at the outer boundary of the environment was chosen to be random-periodic for molecule motion, i.e., a molecule leaving the system at one point was entering the system at another random position of this boundary, where the newly added molecule was given an entirely new identity. At the inner boundary of the cell surface, reflecting boundary conditions were imposed. By applying these realistic boundary conditions, we ensure that the number of molecules in the system is constant during the simulation time.

For a highly realistic implementation of RL binding dynamics, a continuous space representation was used and combined with the neighbor-list method (25, 26) to speed up the detection of interaction partners in this off-lattice approach. Molecules in motion may approach each other and become overlapping. We implemented a push-back procedure, such that the overlap by the moving molecule was reduced to a point contact with the other molecule. Thus, we imposed the condition that molecules cannot penetrate each other and this choice impacts on the effective reaction volume between the molecules.

For reasons of comparison between the different ABM variants, we use the same time step Δt in each simulation, such that changes in the simulation results can be clearly attributed to differences in the receptor morphology, the dimensionality of motion and/or binding valency. To this end, we determine the time step

$\mathrm{\Delta }t=\text{min}\left(\text{min}\left({k_{\mathit{off}}^{micro}}^{-1},{k_{on}^{micro}}^{-1}\right),\text{min}\left(\mathrm{\Delta }{t}_R,\mathrm{\Delta }{t}_L\right)\right),$ (1)

from the smallest considered rate of binding $\left(k_{on}^{micro}\right)$ and dissociation $\left(k_{\mathit{off}}^{micro}\right)$ as well as the smallest time step associated with a diffusion step in space that does not exceed the radius of receptors (Δs_R) and ligands (Δs_L). The time steps of receptors (Δt_R) and ligands (Δt_L) are given by

$\mathrm{\Delta }{t}_{R,L}=\frac{\mathrm{\Delta }{s}_{R,L}^2}{2d{D}_{R,L}}.$ (2)

The simulation algorithm for RL binding dynamics is based on random selection dynamics (5). Each molecule is updated per time step with regard to its diffusion and interaction that are performed in random order applying the acceptance-rejection method (27). A flowchart of the algorithm is shown in Figure 3. For the model system under consideration, i.e., a B cell with a number of BCR in the order 10⁵ and an equal amount of ligands, simulation run times would exceed all limits. In fact, it can be estimated that the ratio of the typical simulation time over the simulated real time becomes as large as 10⁹. Therefore, since the size of the time step is determined by the accurate resolution of molecular motion and interaction, we down-scale the number of molecules and decrease the system size while keeping the molecular concentration constant. The details of the down-scaling procedure are described in Supplementary Material and the associated values are summarized in Table S2 in Supplementary Material. All simulations were performed after down-scaling the number of molecules by a factor s = 10⁻², i.e., reducing the B cell size by a factor 10 and the number of BCR to the order 10³.

FIGURE 3

Figure 3. Flow chart of the ABM simulation algorithm for receptor–ligand binding kinetics. The gray boxes represent operations and are connected by directive arrows depicting the sequence of the ABM simulation algorithm. In each time step Δt, all agents perform diffusive motion and undergo interactions with other agents in random order until simulation time t_∞is reached. Simulations of all ABM variants are shown in Videos S1–S5 in Supplementary Material.

The ABM framework was implemented in the object-oriented programming language C + +.

2.2. Macroscopic Modeling of Receptor–Ligand Binding

Modeling RL binding from a macroscopic point of view can be done in a straightforward fashion using ordinary differential equations (ODE). This approach is appropriate to describe chemical processes where reaction partners occur in large amounts and are homogeneously distributed in the spatial environment. Consequently, ODE models represent time-dependent changes of molecule concentrations in a continuous and deterministic fashion. We considered the binding of receptors (R) and ligands (L) to form a molecular complex (C) as well as their unbinding:

$R+L\underset{k_{\mathit{off}}^{\mathit{macro}}}{\stackrel{k_{on}^{macro}}{\rightleftarrows }}C.$ (3)

Here, $k_{on}^{macro}$ is the reaction rate for binding, $k_{\mathit{off}}^{macro}$ is the dissociation rate and the corresponding association constant K_a is defined by their ratio: $K_a=k_{on}^{macro}∕k_{\mathit{off}}^{macro}$.

The reaction equation (3) was then translated into the coupled system of ODE:

$\frac{dR}{dt}=-k_{on}^{macro}RL+k_{\mathit{off}}^{macro}C,$ (4)

$\frac{dL}{dt}=-k_{on}^{macro}RL+k_{\mathit{off}}^{macro}C,$ (5)

$\frac{dC}{dt}=+k_{on}^{macro}RL-k_{\mathit{off}}^{macro}C.$ (6)

Assuming that initially no molecular complexes exist, C (t = 0) = 0, it follows from the relations R(t) = R(0) – C(t) and L(t) = L(0) – C(t) that it is sufficient to solve the non-linear equation for C(t):

$\frac{dC}{dt}=\alpha C^2-\beta C+\gamma ,$ (7)

where we defined the constants

$\alpha =k_{on}^{macro},$ (8)

$\beta =k_{\mathit{off}}^{macro}+k_{on}^{macro}\left[R\left(0\right)+L\left(0\right)\right],$ (9)

$\gamma =k_{on}^{macro}R\left(0\right)L\left(0\right).$ (10)

The ODE for C(t) can be solved by the separation of variables and yields the analytical solution:

$C\left(t\right)=C_{-}{C}_{+}\frac{1-e^{\alpha \left(C_{+}-C_{-}\right)t}}{C_{-}-C_{+}{e}^{\alpha \left(C_{+}-C_{-}\right)t}}$ (11)

with

$C_{\pm }=\frac{\beta }{2\alpha }\pm \sqrt{\frac{{\beta }^2}{4{\alpha }^2}-\frac{\gamma }{\alpha }}.$ (12)

Note that the concentration C(t) is associated with the number of receptor–ligand (RL) complexes in the microscopic model (see Materials and Methods section 2.1.2).

2.3. Mapping Microscopic and Macroscopic Binding Kinetics

A relation between the macroscopic and microscopic viewpoint on the binding kinetics of receptors and ligands can be established via the corresponding reaction rates for RL binding $k_{on}^{macro}$ and $k_{on}^{micro}$. Given the concentration of molecular complexes C(t) (see equation (11)), we fit this analytical solution from macroscopic binding kinetics to the numerical results of simulations obtained from ABM at the microscopic level. This yields the desired relation $k_{on}^{macro}\left(k_{on}^{micro}\right)$ that can be compared for different ABM variants.

The fitting procedure was performed within the open source programming language R (28). We used the function nls() that returns optimal parameter values of non-linear model equations by least-squares fitting. In particular, we used the fitting algorithm option “port” that refers to the adaptive non-linear least-squares algorithm NL2SOL (29) provided by the Port library. The algorithm adaptively switches between the Gauss-Newton method and an augmented Hessian approximation (30).

In practice, we applied the fitting procedure in two different respects: (i) The macroscopic binding rate $k_{on}^{macro}$ in equation (11) was estimated from fitting to the data points obtained from numerical simulations with the ABM over time. (ii) The values determined for $k_{on}^{macro}$ were used as data points to fit the optimal parameter values of the Hill equation $k_{on}^{macro}\left(k_{on}^{micro}\right)$ (see equation (13)) in order to map the microscopic and macroscopic binding kinetics.

3. Results

In this section, we present our simulation results on receptor–ligand (RL) binding by comparing the dynamics of individual receptors and ligands at the microscopic level with the population kinetics at the macroscopic level. The population kinetics can be straightforwardly described by a coupled system of ordinary differential equations (ODE), whereas agent-based models (ABM) resolve spatial structures of receptors and ligands and account for the dimensionality of the spatial environment in which these molecules diffuse and interact. In particular, we study monovalent receptors with different morphologies, i.e., being either spherically shaped (O) or Y-shaped (Y), and in settings with different dimensionality of motion, i.e., in solution (SOL) or membrane anchored (MEM). While ligands are throughout considered as being in solution and as having spherical shape, the four combinations of receptor properties give rise to four different ABM variants that are denoted by their receptor properties, respectively, as O-SOL, O-MEM, Y-SOL, and Y-MEM. These are schematically depicted in Figure 1 and the differences between receptors are shown in Figure 2. In addition, videos of simulations for the different ABM variants with monovalent receptors are provided in Videos S1–S5 in Supplementary Material, where Videos S1-S4 represent down-scaled systems with factor s = 10⁻², while Video S5 shows a simulation of ABM variant Y-MEM with s = 1. A flow chart of the simulation algorithm is provided in Figure 3 and details on the implementation of the ABM and on the model parameters are given in the Materials and Methods section.

3.1. Binding Kinetics for Different Receptor Properties Qualitatively Comparable

The binding kinetics at the macroscopic level, which can be determined from the analytical solution of the ODE model (see Materials and Methods section), was observed to be in qualitative agreement with the simulation results of all four ABM variants with monovalent receptors at the microscopic level. This can be seen from the ABM simulation results in Figure 4, where the microscopic rate for RL dissociation was fixed at $k_{\mathit{off}}^{micro}=0.1{\text{s}}^{-1}$, while the microscopic rate for RL binding was set to $k_{on}^{micro}={10}^6{\text{s}}^{-1}$ (Figure 4A) and $k_{on}^{micro}={10}^7{\text{s}}^{-1}$ (Figure 4B). Note that we provide the concentration of molecular complexes in units 1/µ m³ to enable the comparison of the binding dynamics simulated by ODE and ABM variants with soluble and membrane-anchored receptors. Since the initial numbers of receptors and ligands as well as the system volumes are identical in all models and simulations, we basically perform a comparison with regard to the number of complexes in each system. In general, we observed that the impact of the stochasticity on RL binding dynamics in the ABM is small, e.g., the relative standard deviation in the number of RL complexes was found to be around 1% for equilibrated systems (see the thickness of curves in pale colors in Figure 4). This is due to the large number of molecules in each simulation, such that five repetitions—involving in total the simulation of 10⁴ molecules—yielded vanishingly small standard deviations.

FIGURE 4

Figure 4. Receptor–ligand binding kinetics for the four ABM variants. Time-dependent concentration of RL complexes with monovalent receptors as obtained from simulations of all four ABM variants with dissociation rate $k_{\mathit{off}}^{micro}=0.1{\text{s}}^{-1}$ and binding rate (A) $k_{on}^{micro}={10}^6{\text{s}}^{-1}$ or (B) $k_{on}^{micro}={10}^7{\text{s}}^{-1}$. Dark and pale lines in different colors represent, respectively, mean values and standard deviations of five simulation runs per ABM variant. Dashed lines indicate the corresponding ODE models after fitting the macroscopic reaction rate $k_{on}^{macro}$.

We generally found a decrease in the concentration of free receptors and ligands with time, which was naturally associated with an increase in the concentration of RL complexes. This observation was robust against variations in the receptor properties, i.e., all four ABM variants—O-SOL, O-MEM, Y-SOL, and Y-MEM—showed the same qualitative behavior. Thus, the qualitative agreement with the macroscopic binding kinetics based on the ODE was not limited to the ABM variant O-SOL as its direct microscopic counterpart. Therefore, in what follows, the analytical ODE solution can be used to fit the simulation results of all four ABM variants and to characterize them by their quantitative differences in the macroscopic binding rate $k_{on}^{macro}$. Note that this is the only free model parameter, since the dissociation of RL complexes occurs spontaneously at both the microscopic and macroscopic level implying that the corresponding rates are identical: $k_{\mathit{off}}^{macro}=k_{\mathit{off}}^{micro}$. Arguments for this relation between macroscopic and microscopic dissociation rates are provided based on the analysis in Supplementary Material.

3.2. Receptor Properties Have Quantitative Impact on Binding Kinetics

At the quantitative level, we observed differences in the binding kinetics depending on the receptor properties as well as on the microscopic binding rate $k_{on}^{micro}$. As could be expected, formation of RL complexes occurred slower for smaller $k_{on}^{micro}={10}^6{\text{s}}^{-1}$ (Figure 4A) than for larger $k_{on}^{micro}={10}^7{\text{s}}^{-1}$ (Figure 4B). Moreover, for a fixed value $k_{on}^{micro}$, the ABM variants with monovalent receptors in solution—O-SOL (red lines) and Y-SOL (blue lines)—exhibited quantitative agreement in the binding kinetics. While for the corresponding ABM variants with membrane-anchored receptors—O-MEM (orange lines) and Y-MEM (green lines)—this quantitative agreement was also observed, a quantitative difference in the binding kinetics between receptors in solution and membrane-anchored receptors was clearly visible (see Figure 4).

Using the analytical ODE solution of the binding kinetics, we fitted the simulation results of all four ABM variants to characterize them by their quantitative differences in the macroscopic binding rate $k_{on}^{macro}$. The fitted curves are shown in Figure 4 and yielded for $k_{on}^{micro}={10}^6{\text{s}}^{-1}$ (Figure 4A) the values $k_{on}^{macro}\approx 1.9\mathrm{\mu }{\text{m}}^3{\text{s}}^{-1}$ for the ABM variants O-SOL and Y-SOL and $k_{on}^{macro}\approx 0.6\mathrm{\mu }{\text{m}}^3{\text{s}}^{-1}$ for the ABM variants O-MEM and Y-MEM. For $k_{on}^{micro}={10}^7{\text{s}}^{-1}$ (Figure 4B), we obtained the values $k_{on}^{macro}\approx 10.5\mathrm{\mu }{\text{m}}^3{\text{s}}^{-1}$ for the ABM variants O-SOL and Y-SOL and $k_{on}^{macro}\approx 1.7\mathrm{\mu }{\text{m}}^3{\text{s}}^{-1}$ for the ABM variants O-MEM and Y-MEM. It should be noted that the goodness of the fit, which was evaluated by the error of least squares fitting, was comparable for all simulations with microscopic binding rates in the range ${10}^4{\text{s}}^{-1}\le k_{on}^{micro}\le {10}^6{\text{s}}^{-1}$. Even though for $k_{on}^{micro}>{10}^6$ the error of least squares fitting for ABM variants with membrane-anchored receptors can be up to two orders of magnitude larger than for those with receptors in solution (see Figure S1 in Supplementary Material), all fitted curves still represented a fair representation of the simulation results (see Figure 4B).

These results were the first indication that the receptor morphology plays a relatively minor role in the binding kinetics compared to the dimensionality of motion of receptors, i.e., whether receptors diffuse in three-dimensional solution or on the surface of a cell. To further analyze these findings, we decided to establish a detailed quantitative mapping between the macroscopic and microscopic binding rates.

3.3. Quantitative Mapping of the Macroscopic and Microscopic Binding Rates Reveals Impact of Dimensionality of Motion

We performed numerical simulations to quantify the difference in monovalent RL binding as a function of receptor properties. All four ABM variants were applied using the fixed dissociation rate $k_{\mathit{off}}^{micro}=k_{\mathit{off}}^{macro}=0.1{\text{s}}^{-1}$ and varying the microscopic binding rate in the range ${10}^4{\text{s}}^{-1}\le k_{on}^{micro}\le 2.5\times {10}^7{\text{s}}^{-1}$. The corresponding macroscopic binding rate $k_{on}^{macro}$ was determined for each numerical experiment from the best fit of the analytical solution of the ODE model to the simulation result of the ABM. The resulting function $k_{on}^{macro}\left(k_{on}^{micro}\right)$ is shown in Figure 5 for each ABM variant. The steady state concentrations of complexes and receptors obtained by fitting the ODE kinetics to the dynamics of the four various ABM variants are summarized in Tables S3–S6 in Supplementary Material.

FIGURE 5

Figure 5. Mapping of microscopic and macroscopic binding rates for different ABM variants. Simulation of all four ABM variants for varying $k_{on}^{micro}$ and the fitted reaction rate $k_{on}^{\mathit{macro}}$ of the ODE models. Solid lines represent Hill functions with parameters fitted to the data points $k_{on}^{\mathit{macro}}\left(k_{on}^{micro}\right)$. Results for ABM variants are similar for the same dimensionality of motion for receptors, i.e., either in solution (O-SOL, Y-SOL) or membrane anchored (O-MEM, Y-MEM), but are distinct for ABM variants with soluble and membrane-anchored receptors. The dotted line represents the binding rate as determined by Collins and Kimball (see equation (14)) that is, as expected, comparable to the simulation result for ABM variant O-SOL.

As expected from our previous considerations, the quantitative difference between morphologies of monovalent receptors is negligible compared to the dimensionality of motion, i.e., whether receptors were diffusing in solution or within the membrane on the surface of a cell. Moreover, the numerical results $k_{on}^{macro}\left(k_{on}^{micro}\right)$ in Figure 5 resemble Hill functions,

$k_{on}^{macro}\left(k_{on}^{micro}\right)=\frac{a{k}_{on}^{micro}}{b+k_{on}^{micro}},$ (13)

with parameters a and b that are specific for given receptor properties. Here, a denotes the upper limit for the macroscopic binding rate, $k_{on}^{macro}\left(k_{on}^{micro}\gg b\right)\to a$, and b is a constant that determines the slope of the Hill function, $k_{on}^{macro}\left(k_{on}^{micro}\ll b\right)\to \left(a∕b\right)k_{on}^{micro}$, while at intermediate value $k_{on}^{micro}=b$ the Hill function attains half of its maximal value: $k_{on}^{macro}\left(k_{on}^{micro}=b\right)=a∕2$. The two parameters can be determined from a fit to the numerical simulations and the resulting curves are shown in Figure 5 as solid lines. The corresponding values are summarized in Table S7 in Supplementary Material for the four ABM variants.

The observed functional dependence of $k_{on}^{macro}$ on $k_{on}^{micro}$ is in agreement with theoretical considerations by Collins and Kimball on binding reactions of diffusing receptors and ligands in three spatial dimensions (31–33). They arrived at the expression

$k_{on}^{macro}\left(\mathrm{\kappa }\right)=\frac{k_s\mathrm{\kappa }}{k_s+\mathrm{\kappa }},$ (14)

where $k_s=4\pi \left(r_L+r_R\right)\left(D_L+D_R\right)$ denotes the diffusion-controlled reaction rate that was previously introduced by Von Smoluchowski (34) and that depends on the radii of receptor (r_R) and ligand (r_L) as well as on the diffusion coefficients of receptor (D_R) and ligand (D_L). This rate refers to the frequency at which diffusing receptors and ligands come into contact, i.e., have the distance r_R + r_L. Furthermore, κ denotes the intrinsic reaction rate, $\mathrm{\kappa }=V_r{k}_{on}^{micro}$, which is directly related to the microscopic binding rate $k_{on}^{micro}$ and the reaction volume $V_r=\left(4∕3\right)\mathrm{\pi }{\left(r_L+r_R\right)}^3$ (35, 36). Combining equations (13) and (14) yields the following relationships:

$a=k_s,$ (15)

$b=\frac{k_s}{V_r}.$ (16)

It should be stressed that this correspondence can strictly speaking only be applied to monovalent receptors with spherical morphology and to RL binding in three-dimensional solution with receptor and ligand being allowed to penetrate each other. In other words, equations (15) and (16) could only be expected to hold for the ABM variant O-SOL, however, even this scenario is different from the theoretical considerations in that molecules are not allowed to penetrate each other in our ABM. In the ABM, we generally do not allow for molecular penetration in RL interactions, which reduces their possible overlap to a point contact. The implementation of push-back collisions between molecules effectively reduces the reaction volume V_r, i.e., we set V_r → f_rV_r with scaling factor f_r ≤ 1. This parameter will only affect the slope of the Hill function, while it was observed in Figure 5 that the upper limit of the macroscopic binding rate, k_s, does as well depend on the receptor properties. To account for these observations, we set k_s → f_sk_s with scaling factor f_s. It then follows that f_r and f_s can be computed from the equations

$f_s=\frac{a}{k_s},$ (17)

$f_r=\frac{f_s{k}_s}{b{V}_r}$ (18)

in terms of the two fitting parameters a and b (see Table S7 in Supplementary Material). The resulting scaling factors are summarized in Table S8 in Supplementary Material.

As could be expected, for the ABM variant O-SOL we found the scaling factor $f_s^{O-SOL}=1.02$ to be close to 1, implying that the upper limit for the macroscopic binding rate as predicted by Collins and Kimball was quantitatively recovered (31–33). Regarding the increase of $k_{on}^{macro}$ as a function of $k_{on}^{micro}$, we found the difference in the underlying assumptions on RL interactions to be reflected by a decrease in the reaction volume V_r with scaling factor $f_r^{O-SOL}=0.79$.

We compared the scaling factors for the other ABM variants and present the results relative to ABM variant O-SOL in Figure 6. The scaling factor $f_r^{Y-SOL}$ of ABM variant Y-SOL was found to be similar to $f_r^{O-SOL}$ with a relative decrease of only 4%, whereas this scaling factor for the ABM variants with membrane-anchored receptors, i.e., $f_r^{O-MEM}$ and $f_r^{Y-MEM}$, was decreased by 74 and 61%, respectively. Furthermore, as shown in Figure 6, the scaling factors $f_s^{O-MEM}$ and $f_s^{Y-MEM}$ for membrane-anchored receptors were found to be decreased from $f_s^{O-SOL}$ by 77 and 69%, respectively, indicating a significant change in the upper limit of the macroscopic binding rate. On the other hand, this scaling factor was always somewhat higher for membrane-anchored receptors, i.e., ABM variants O-MEM and Y-MEM, compared to their respective counterparts with soluble receptors.

FIGURE 6

Figure 6. Quantitative difference in the scaling factors of ABM variants relative to O-SOL. The scaling factors f_s and f_r are calculated from equations (17) and (18) for parameters specific to the considered ABM variant. ABM variant O-SOL resembles the conditions of the theoretical considerations by Collins and Kimball (31–33) most of all. Scaling factors of ABM variants with membrane-anchored receptors that are either spherically shaped (O-MEM) or Y-shaped (Y-MEM) exhibit similar but clear differences to ABM variant O-SOL, whereas ABM variant Y-SOL is most similar to O-SOL.

We checked the dependency of the mapping between macroscopic and microscopic binding rates (see Figure 5) as well as the scaling factors f_s and f_r (see Figure 6) on the down-scaling factor s of the simulated ABM variants. It was generally observed that simulations for soluble receptors were not affected by the system down-scaling, whereas in simulations for membrane-anchored receptors increasing the down-scaling factor s resulted into lower values for $k_{on}^{macro}$ as a function of $k_{on}^{micro}$. This implies that the difference between ABM variants with soluble and membrane-anchored receptors as observed in Figure 5 as well as the distances between the respective scaling factors in Figure 6 represents a lower limit.

Since the diffusion coefficients of receptors in the soluble (D_R = 90 µm² s⁻¹) and membrane-anchored (D_R = 0.05 µm² s⁻¹) variant differed by orders of magnitude, we checked whether differences in the upper limit of the macroscopic binding rate were indeed merely a consequence of the dimensionality of motion rather than of the magnitude of the diffusion coefficient itself. This was done by running simulations with interchanged diffusion coefficients, i.e., ABM variant O-SOL with D_R = 0.05 µm² s⁻¹ and ABM variant O-MEM with D_R = 90 µm² s⁻¹. However, even this dramatic modification of diffusion coefficients did not eliminate the significant difference in the dependence of $k_{on}^{macro}$ on $k_{on}^{micro}$ between the ABM variants (see Figures S2 and S3 in Supplementary Material).

Taken together, our quantitative analysis of monovalent RL binding kinetics revealed the impact of receptor properties on the macroscopic binding rate and by that on the association constant of the RL binding. It was shown that the diffusion coefficients of receptors and their morphology have minor effects, whereas the strongest impact was due to the dimensionality of motion. Compared to soluble receptors in three dimensions, RL binding kinetics of membrane-anchored receptors on a cellular surface were retarded and could not achieve comparably high association constants. In what follows, we consider the impact of the binding valency by taking into account that the Y-shaped receptors can bind a ligand at each receptor arm.

3.4. Binding Valency Reduces Differences in the Binding Kinetics of BCR and Antibodies

To investigate the influence of the receptor binding valency on the binding kinetics for monovalent receptors (see Figure 5), we modified ABM variants Y-MEM and Y-SOL as to allow for bivalent binding of the Y-shaped receptors, i.e., a ligand can bind at each of the two receptor arms. Thus, in these ABM variants the term complex refers to receptors that are bound to either one or two ligands. The simulations were performed with varied binding rate $k_{on}^{micro}$ between 5 × 10⁶ and 2.5 × 10⁷ s⁻¹. The temporal course of the binding kinetics for simulations of the bivalent and monovalent ABM variants is shown in Figure 7. The simulations of $k_{on}^{micro}=1\times {10}^7{\text{s}}^{-1}$ exhibit the typical relations between the binding kinetics of the ABM variants. As could be expected, both ABM variants with bivalent receptors showed a faster binding kinetics and also reached higher association constants than their monovalent counterparts. In Figure 8, we show the relative difference in receptor-bound ligands for ABM variant Y-MEM relative to ABM variant Y-SOL and for different values of $k_{on}^{micro}$. This difference is significantly smaller (down to 72%) for bivalent receptors compared to monovalent receptors, and in the limit of long times this difference vanishes only for bivalent but not for monovalent receptors. These results indicate that the binding valency makes a clear difference for RL binding: In the case of monovalent receptors, the dimensionality of motion induces a significant difference in the binding kinetics, whereas this difference is largely compensated by the bivalency of receptors. Thus, it turns out that membrane-anchored BCR and soluble antibodies do reach comparable association constants for bivalent receptors.

FIGURE 7

Figure 7. Kinetics of bound ligands for Y-MEM and Y-SOL ABM variants with either monovalent or bivalent receptors. Time-dependent concentration of bound ligands for ABM variants Y-MEM and Y-SOL for models with either monovalent receptors or bivalent receptors. All models were simulated with dissociation rate $k_{\mathit{off}}^{micro}=0.1{\text{s}}^{-1}$ and binding rate $k_{on}^{micro}={10}^7{\text{s}}^{-1}$. Dark and pale lines in different colors represent, respectively, mean values and standard deviations of five simulation runs per ABM variant.

FIGURE 8

Figure 8. Relative differences between ABM variants Y-MEM and Y-SOL with either monovalent or bivalent receptors. Temporal evolution of the relative differences of bound ligands between ABM variants Y-MEM and Y-SOL for models with either monovalent receptors (monovalent ABM variant) or bivalent receptors (bivalent ABM variant). The colors refer to ABM variants with varying binding rates $k_{on}^{micro}$.

In order to investigate whether these observations are caused by the effectively twofold number of binding sites for the bivalent receptors, we performed simulations with ABM variants that have twice as much monovalent receptors than the so far applied physiological number of receptors $\left(N_p^R\right)$. The binding kinetics of ABM variants with $N^R=2\times N_p^R$ monovalent receptors turned out to be even faster as the binding kinetics of bivalent ABM variants with $N^R=N_p^R$ (see Figure S4 in Supplementary Material). Additionally, the relative differences between binding kinetics of ABM variants with soluble and membrane-bound receptors vanishes with increasing time, and this occurs slightly faster as for ABM variants with bivalent receptors (see Figure S5 in Supplementary Material). These results indicate that comparable association constants of membrane anchored and soluble receptors can be observed for systems with higher amounts of binding sites at receptors.

4. Discussion

The focus of this study on receptor–ligand (RL) binding was twofold. Firstly, we established a quantitative mapping between macroscopic binding rates of an ordinary differential equation (ODE) model and their microscopic equivalents as obtained from simulating the spatiotemporal binding kinetics by agent-based models (ABM). Secondly, we investigated the impact of various properties of B cell-derived receptors—such as their dimensionality of motion, morphology and binding valency—on the RL binding kinetics.

Regarding the quantitative mapping of binding rates, we recovered for fixed dissociation rates $k_{\mathit{off}}^{micro}=k_{\mathit{off}}^{macro}=0.1{\text{s}}^{-1}$ the non-linear relationship between the binding rates $k_{on}^{macro}$ and $k_{on}^{micro}$. This resembles a Hill-type function (see Figure 5), which is in line with theoretical predictions by Collins and Kimball (31–33). Scanning $k_{on}^{micro}$ over more than four orders of magnitude, we obtained upper limiting values for $k_{on}^{macro}$ in the range 10⁰–10¹ µm³ s⁻¹, which corresponds to 10⁸–10⁹ M⁻¹ s⁻¹ using Avogadro’s number. For $k_{\mathit{off}}^{macro}=0.1{\text{s}}^{-1}$, the resulting association constant is K_a = 10¹⁰ M⁻¹. This is in agreement with experimentally measured values for BCR-antigen binding, where typical values up to K_a = 10¹⁰ M⁻¹ are reached (37, 38), which is a strong indication for our ABM variants to be realistic and quantitative to-scale representations of RL binding.

The ABM variants were implemented in three-dimensional representations of continuous space and RL binding was simulated by the random selection method (5). We implemented different ABM variants where binding of spherical ligands occurs either with soluble receptors or with membrane-anchored receptors. The receptors are either spherically shaped or Y-shaped and can be mono- or bivalent. We simulated RL binding in identical environments to allow for quantitative comparisons of the different scenarios. In particular, we considered the Y-shaped and bivalent antibodies in solution and the B cell receptors (BCR) as their membrane-anchored counterparts on a spherical cell to be an appropriate example. In previous work on BCR binding, ABM implementations typically involved simplifications with regard to the spatial representation, i.e., using a planar cell surface and imposing a spatial grid for molecule diffusion (39, 40) and have been applied to simulate the immunological synapse involving B cells (41–45) or T cells (46, 47). Besides this work on immune cell receptor–ligand interaction, there exist software packages for the simulation of various type, such as Smoldyn (48) and MCell (49, 50). Even though these simulators represent molecular diffusion in lattice-free continuous space, they lack features that are essential in the present study. For example, Smoldyn represents molecules in a point-like fashion (48, 51–53), while MCell does only allow to determine an upper limit of the simulation time step Δt (54) implying that simulations with different model systems may differ in the time step Δt. Therefore, we did not consider these simulators suitable for the investigation of morphological aspects of receptors and for comparing models at the microscopic and macroscopic scale. Moreover, the RL binding of soluble and membrane-anchored receptors was previously also investigated by non-spatial ODE models (55, 56). These two-step ODE models comprise the process of encounter formation by molecule diffusion and the reaction process itself, so that molecular parameters, like diffusion constant and size, could also be incorporated. However, several simplifications were made, such as the derivation of the binding rate of membrane-bound receptors from cell–ligand interaction rate, which turned out to be not applicable in general (55, 56).

To study the impact of various receptor properties on RL binding kinetics, we compared scenarios that differ in the dimensionality of motion, morphology and binding valency of receptors. These receptor properties were investigated since they are characteristic for B cell-derived receptors that play a key role in the adaptive immune response. Interestingly, the RL binding kinetics for monovalent Y-shaped receptors was observed to be quantitatively comparable to that of spherical receptors (see Figure 5), i.e., the difference in the morphology of monovalent receptors did not reveal a substantial impact. In contrast, the dimensionality of motion for BCR compared to soluble antibodies did reveal a clear difference in the binding kinetics, i.e., the association constants were found to be significantly lower for membrane-anchored receptors compared to soluble receptors (see Figure 5). Furthermore, our results show that the diffusion constant of receptors, which is much smaller for membrane-anchored molecules as for soluble molecules, does not strongly influence the observed differences in the binding kinetics. This suggest that the difference in the association constants for soluble and membrane-anchored monovalent receptors originate from the difference in the dimensionality of motion. However, this difference was largely compensated by taking into account that BCR and soluble antibodies are bivalent (see Figure 8), i.e., the relative difference in the binding kinetics of membrane-anchored and soluble receptors vanished only in the case of bivalent receptors. It is generally known that the bivalency of BCR supports cross-linking in the binding to multivalent ligands. However, the current findings suggest that the bivalency of BCR does also compensate the difference in the association constant that exist for monovalent receptors between the soluble and membrane-anchored variants.

In the future, the extensibility of the current simulation framework can be exploited to study more complex scenarios. For example, antigens may be represented by multivalent ligands that do not only allow for cross-linking of BCR but also binding to coreceptors required for B cell activation. This enables to study the important process of BCR clustering on the cell surface (57–59) that has also been the subject of theoretical investigations (39, 40, 60, 61). We envisage that such studies will strongly benefit from an image-based systems biology approach, for example, as applied by Mech et al. (62) and conceptionally reviewed by Medyukhina et al. (63). Recently, we took the first steps toward an image-based investigation of B cell activation that requires the concerted action of various receptors and ligands (64). Based on these data, our ABM can be extended by various agent types with specific properties to predict prerequisites for experimentally observed molecular patterns. Moreover, the ABM variants could be modified to represent various receptor properties of different antibody isotypes and/or subclasses, which would allow investigating the impact of specific receptor properties on the RL binding kinetics. Based on this modification, the impact of naturally occurring antibody complexes, such as IgA dimers and IgM pentamers, could be investigated. Furthermore, extending the ABM to represent arbitrarily shaped cells that are brought in close contact, it can be used to simulate the molecular patterns during synapse formation involving B cells, T cells as well as phagocytes (65–68). This would enable to investigate the impact of the dimensionality of motion of ligands that is reported to be an important parameter for regulating B cell activation and signaling (69).

Author Contributions

TL and MF conceived and designed the study, evaluated and analyzed the results, and wrote the manuscript and critically revised it. MF contributed materials and computational resources. TL processed the data, implemented, and applied the computational algorithm.

Conflict of Interest Statement

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

Acknowledgments

We acknowledge support by Bertram Vogel regarding the initial implementation of the agent-based model.

Funding

This work was financially supported by the Center for Sepsis Control and Care (CSCC) (Project Quantim to MTF, FKZ 01EO1502) that is funded by the Federal Ministry for Education and Research (BMBF) and by the CRC/TR124 FungiNet (Project B4 to MTF) that is funded by the Deutsche Forschungsgemeinschaft (DFG).

Supplementary Material

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

Video S1. Simulation of down-scaled scaled O-SOL ABM variant. ABM simulation with monovalent receptors (blue objects) that are spherically shaped and move in solution by performing three-dimensional diffusion. Upon contact between receptors and ligands (red objects) these may bind and form RL complexes (green objects) depending on the binding rate $k_{on}^{micro}$ = 2.5 × 10⁷ s⁻¹. The system is down-scaled with factor s = 0.01 (see Supplementary Material) and values of model parameters are provided in Tables S1 and S2 in Supplementary Material. The video is composed of 15 frames s⁻¹ and the simulation time between two consecutive frames is 6.8 × 10⁻⁸ s. A high-resolution video is available for download from https://asbdata.hki-jena.de/LehnertFigge2017_FrontImmun/.

Video S2. Simulation of down-scaled O-MEM ABM variant. ABM simulation with monovalent receptors (blue objects) that are spherically shaped and move in the cell membrane by performing two-dimensional diffusion. Upon contact between receptors and ligands (red objects) these may bind and form RL-complexes (green objects) depending on the binding rate $k_{on}^{micro}$ = 2.5 × 10⁷ s⁻¹. The system is down-scaled with factor s = 0.01 (see Supplementary Material) and values of model parameters are provided in Tables S1 and S2 in Supplementary Material. The video is composed of 15 frames s⁻¹ and the simulation time between two consecutive frames is 6.8 × 10⁻⁸ s. A high-resolution video is available for download from https://asbdata.hki-jena.de/LehnertFigge2017_FrontImmun/.

Video S3. Simulation of down-scaled Y-SOL ABM variant. ABM simulation with monovalent receptors (blue objects) that are Y-shaped and move in solution by performing three-dimensional diffusion. Upon contact between receptors and ligands (red objects) these may bind and form RL-complexes (green objects) depending on the binding rate $k_{on}^{micro}$ = 2.5 × 10⁷ s⁻¹. The system is down-scaled with factor s = 0.01 (see Supplementary Material) and values of model parameters are provided in Table S1 in Supplementary Material and Supplementary Material. The video is composed of 15 frames s⁻¹ and the simulation time between two consecutive frames is 6.8 × 10⁻⁸ s. A high-resolution video is available for download from https://asbdata.hki-jena.de/LehnertFigge2017_FrontImmun/.

Video S4. Simulation of down-scaled Y-MEM ABM variant. ABM simulation with monovalent receptors (blue objects) that are Y-shaped and move in the cell membrane by performing two-dimensional diffusion. Upon contact between receptors and ligands (red objects) these may bind and form RL-complexes (green objects) depending on the binding rate $k_{on}^{micro}$ = 2.5 × 10⁷ s⁻¹. The system is down-scaled with factor s = 0.01 (see Supplementary Material) and values of model parameters are provided in Tables S1 and S2 in Supplementary Material. The video is composed of 15 frames s⁻¹ and the simulation time between two consecutive frames is 6.8 × 10⁻⁸ s. A high-resolution video is available for download from https://asbdata.hki-jena.de/LehnertFigge2017_FrontImmun/.

Video S5. Simulation of Y-MEM ABM variant. ABM simulation with monovalent receptors (blue objects) that are Y-shaped and move in the cell membrane by performing two-dimensional diffusion. Upon contact between receptors and ligands (red objects) these may bind and form RL-complexes (green objects) depending on the binding rate $k_{on}^{micro}$ = 2.5 × 10⁷ s⁻¹. The system is down-scaled with factor s = 0.01 (see Supplementary Material) and values of model parameters are provided in Table S1 in Supplementary Material. The video is composed of 15 frames s⁻¹ and the simulation time between two consecutive frames is 6.8 × 10⁻⁸ s. A high-resolution video is available for download from https://asbdata.hki-jena.de/LehnertFigge2017_FrontImmun/.

References

1. Resat H, Petzold L, Pettigrew MF. Kinetic modeling of biological systems. Methods Mol Biol (2009) 541:311–35. doi:10.1007/978-1-59745-243-4_14

CrossRef Full Text | Google Scholar

2. Faro J, Castro M, Molina-París C. A unifying mathematical framework for experimental TCR-pMHC kinetic constants. Sci Rep (2017) 7:46741. doi:10.1038/srep46741

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Andrews SS, Arkin AP. Simulating cell biology. Curr Biol (2006) 16(14):R523–7. doi:10.1016/j.cub.2006.06.048

CrossRef Full Text | Google Scholar

4. Goldstein B, Faeder JR, Hlavacek WS. Mathematical and computational models of immune-receptor signalling. Nat Rev Immunol (2004) 4(6):445–56. doi:10.1038/nri1374

CrossRef Full Text | Google Scholar

5. Figge MT. Stochastic discrete event simulation of germinal center reactions. Phys Rev E Stat Nonlin Soft Matter Phys (2005) 71(5):051907. doi:10.1103/PhysRevE.71.051907

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J Phys Chem (1977) 81(25):2340–61. doi:10.1021/j100540a008

CrossRef Full Text | Google Scholar

7. Gillespie DT. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys (1976) 22:403–34. doi:10.1016/0021-9991(76)90041-3

CrossRef Full Text | Google Scholar

8. Gibson MA, Bruck J. Efficient exact stochastic simulation of chemical systems with many species and many channels. J Phys Chem A (2000) 104(9):1876–89. doi:10.1021/jp993732q

CrossRef Full Text | Google Scholar

9. Yu JS, Bagheri N. Multi-class and multi-scale models of complex biological phenomena. Curr Opin Biotechnol (2016) 39:167–73. doi:10.1016/j.copbio.2016.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Bonabeau E. Agent-based methods and techniques for simulating human systems. Proc Natl Acad Sci U S A (2002) 99(10):7280–7. doi:10.1073/pnas.082080899

CrossRef Full Text | Google Scholar

11. Takahashi K, Vel Arjunan SN, Tomita M. Space in systems biology of signaling pathways – towards intracellular molecular crowding in silico. FEBS Lett (2005) 579(8):1783–8. doi:10.1016/j.febslet.2005.01.072

CrossRef Full Text | Google Scholar

12. Liu W, Meckel T, Tolar P, Sohn HW, Pierce SK. Antigen affinity discrimination is an intrinsic function of the B cell receptor. J Exp Med (2010) 207(5):1095–111. doi:10.1084/jem.20092123

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Coico R, Sunshine G. Immunology: A Short Course. 7th ed. Oxford, UK: John Wiley & Sons Ltd (2015).

Google Scholar

14. Figge MT, Garin A, Gunzer M, Kosco-Vilbois M, Toellner K-M, Meyer-Hermann M. Deriving a germinal center lymphocyte migration model from two-photon data. J Exp Med (2008) 205(13):3019–29. doi:10.1084/jem.20081160

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Meyer-Hermann M, Figge MT, Toellner KM. Germinal centres seen through the mathematical eye: B-cell models on the catwalk. Trends Immunol (2009) 30(4):157–64. doi:10.1016/j.it.2009.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Garin A, Meyer-Hermann M, Contie M, Figge MT, Buatois V, Gunzer M, et al. Toll-like receptor 4 signaling by follicular dendritic cells is pivotal for germinal center onset and affinity maturation. Immunity (2010) 33(1):84–95. doi:10.1016/j.immuni.2010.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Raychaudhuri S. The problem of antigen affinity discrimination in B-cell immunology. ISRN Biomath (2013) 2013:1–18. doi:10.1155/2013/845918

CrossRef Full Text | Google Scholar

18. Zhang Y, Meyer-Hermann M, George LA, Figge MT, Khan M, Goodall M, et al. Germinal center B cells govern their own fate via antibody feedback. J Exp Med (2013) 210(3):457–64. doi:10.1084/jem.20120150

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Horn F, Heinekamp T, Kniemeyer O, Pollmächer J, Valiante V, Brakhage AA. Systems biology of fungal infection. Front Microbiol (2012) 3:108. doi:10.3389/fmicb.2012.00108

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zhang C, Chen J, DeLisi C. Protein-protein recognition: exploring the energy funnels near the binding sites. Proteins (1999) 34(2):255–67. doi:10.1002/(SICI)1097-0134(19990201)34:2<255::AID-PROT10>3.0.CO;2-O

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Tsai C-J, Kumar S, Ma B, Nussinov R. Folding funnels, binding funnels, and protein function. Protein Sci (1999) 8(6):1181–90. doi:10.1110/ps.8.6.1181

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Northrup SH, Erickson HP. Kinetics of protein-protein association explained by Brownian dynamics computer simulation. Proc Natl Acad Sci U S A (1992) 89:3338–42.

Google Scholar

23. Tovchigrechko A, Vakser IA. How common is the funnel-like energy landscape in protein–protein interactions? Protein Sci (2001) 10(8):1572–83. doi:10.1110/ps.8701

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Einstein A. Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Macromol Symp (1905) 322(8):549–60.

Google Scholar

25. Lehnert T, Timme S, Pollmächer J, Hünniger K, Kurzai O, Figge MT. Bottom-up modeling approach for the quantitative estimation of parameters in pathogen-host interactions. Front Microbiol (2015) 6:608. doi:10.3389/fmicb.2015.00608

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Rapaport DC. The Art of Molecular Dynamics Simulation. New York, NY, USA: Cambridge University Press (2004).

Google Scholar

27. Press W, Teukolsky S, Vetterling W, Flannery B, Ziegel E, Press W, et al. Numerical Recipes: The Art of Scientific Computing. 3rd ed. New York: Cambridge University Press (2007).

Google Scholar

28. Ihaka R, Gentleman R. R: A Language for Data Analysis and Graphics. J Comput Graph Stat (1996) 5(3):299–14. doi:10.2307/1390807

CrossRef Full Text | Google Scholar

29. Dennis JE, Gay DM, Welsch RE. Algorithm 573: NL2SOL – an adaptive nonlinear least-squares algorithm [E4]. ACM Trans Math Software (1981) 7(3):369–83. doi:10.1145/355958.355966

CrossRef Full Text | Google Scholar

30. Butcher J, Jackiewicz Z, Mittelmann H. A nonlinear optimization approach to the construction of general linear methods of high order. J Comput Appl Math (1997) 81(97):181–96. doi:10.1016/S0377-0427(97)00039-3

CrossRef Full Text | Google Scholar

31. Collins FC, Kimball GE. Diffusion-controlled reaction rates. J Colloid Sci (1949) 4(4):425–37. doi:10.1016/0095-8522(49)90023-9

CrossRef Full Text | Google Scholar

32. Collins FC, Kimball GE. Diffusion-controlled reactions in liquid solutions. Indus Eng Chem (1949) 41:2551–3. doi:10.1021/ie50479a040

CrossRef Full Text | Google Scholar

33. Shoup D, Szabo A. Role of diffusion in ligand binding to macromolecules and cell-bound receptors. Biophys J (1982) 40(1):33–9. doi:10.1016/S0006-3495(82)84455-X

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Von Smoluchowski M. Versuch einer mathematischen Theorie der Koagulationskinetik. Phys Chem (1917) 92(1912):156.

Google Scholar

35. Klann MT, Lapin A, Reuss M. Agent-based simulation of reactions in the crowded and structured intracellular environment: influence of mobility and location of the reactants. BMC Syst Biol (2011) 5(1):71. doi:10.1186/1752-0509-5-71

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Klann MT, Koeppl H. Spatial simulations in systems biology: from molecules to cells. Int J Mol Sci (2012) 13(12):7798–827. doi:10.3390/ijms13067798

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Batista FD, Neuberger MS. Affinity dependence of the B cell response to antigen: a threshold, a ceiling, and the importance of off-rate. Immunity (1998) 8(6):751–9. doi:10.1016/S1074-7613(00)80580-4

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Carrasco YR, Fleire SJ, Cameron T, Dustin ML, Batista FD. LFA-1/ICAM-1 interaction lowers the threshold of B cell activation by facilitating B cell adhesion and synapse formation. Immunity (2004) 20(5):589–99. doi:10.1016/S1074-7613(04)00105-0

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Reddy S, Tsourkas PK, Raychaudhuri S. Monte Carlo study of B-cell receptor clustering mediated by antigen crosslinking and directed transport. Cell Mol Immunol (2011) 8(3):255–64. doi:10.1038/cmi.2011.3

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Reddy S, Chilukuri S, Raychaudhuri S. The network of receptors characterize B cell receptor micro- and macroclustering in a Monte Carlo model. J Phys Chem B (2010) 114(1):487–94. doi:10.1021/jp9079074

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Tsourkas PK, Baumgarth N, Simon SI, Raychaudhuri S. Mechanisms of B-cell synapse formation predicted by Monte Carlo simulation. Biophys J (2007) 92(12):4196–208. doi:10.1529/biophysj.106.094995

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Tsourkas PK, Longo ML, Raychaudhuri S. Monte Carlo study of single molecule diffusion can elucidate the mechanism of B cell synapse formation. Biophys J (2008) 95(3):1118–25. doi:10.1529/biophysj.107.122564

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Tsourkas PK, Raychaudhuri S. Modeling of B cell synapse formation by Monte Carlo simulation shows that directed transport of receptor molecules is a potential formation mechanism. Cell Mol Bioeng (2010) 3(3):256–68. doi:10.1007/s12195-010-0123-1

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Tsourkas PK, Liu W, Das SC, Pierce SK, Raychaudhuri S. Discrimination of membrane antigen affinity by B cells requires dominance of kinetic proofreading over serial engagement. Cell Mol Immunol (2012) 9(1):62–74. doi:10.1038/cmi.2011.29

CrossRef Full Text | Google Scholar

45. Tsourkas PK, Somkanya CD, Yu-Yang P, Liu W, Pierce SK, Raychaudhuri S. Formation of BCR oligomers provides a mechanism for B cell affinity discrimination. J Theor Biol (2012) 307:174–82. doi:10.1016/j.jtbi.2012.05.008

CrossRef Full Text | Google Scholar

46. Figge MT, Meyer-Hermann M. Geometrically repatterned immunological synapses uncover formation mechanisms. PLoS Comput Biol (2006) 2(11):e171. doi:10.1371/journal.pcbi.0020171

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Figge MT, Meyer-Hermann M. Modeling receptor-ligand binding kinetics in immunological synapse formation. Eur Phys J D (2009) 51(1):153–60. doi:10.1140/epjd/e2008-00087-1

CrossRef Full Text | Google Scholar

48. Andrews SS, Bray D. Stochastic simulation of chemical reactions with spatial resolution and single molecule detail. Phys Biol (2004) 1(3–4):137–51. doi:10.1088/1478-3967/1/3/001

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Stiles JR, Van Helden D, Bartol TM, Salpeter EE, Salpeter MM. Miniature endplate current rise times less than 100 microseconds from improved dual recordings can be modeled with passive acetylcholine diffusion from a synaptic vesicle. Proc Natl Acad Sci U S A (1996) 93(12):5747–52. doi:10.1073/pnas.93.12.5747

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Kerr RA, Bartol TM, Kaminsky B, Dittrich M, Chang J-CJ, Baden SB, et al. Fast Monte Carlo simulation methods for biological reaction-diffusion systems in solution and on surfaces. SIAM J Sci Comput (2008) 30(6):3126. doi:10.1137/070692017

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Andrews SS. Spatial and stochastic cellular modeling with the Smoldyn simulator. Methods Mol Biol (2012) 804(1):519–42. doi:10.1007/978-1-61779-361-5_26

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Andrews SS. Serial rebinding of ligands to clustered receptors as exemplified by bacterial chemotaxis. Phys Biol (2005) 2(2):111–22. doi:10.1088/1478-3975/2/2/004

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Andrews SS, Addy NJ, Brent R, Arkin AP. Detailed simulations of cell biology with Smoldyn 2.1. PLoS Comput Biol (2010) 6(3):e1000705. doi:10.1371/journal.pcbi.1000705

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Burrage K, Burrage PM, Leier A, Marquez-Lago T, Nicolau DV Jr. Stochastic simulation of spatial modelling of dynamic processes in a living cell. In: Koeppl H, Densmore D, Setti M, Di Bernardo M, editor. Design and Analysis of Biomolecular Circuits. New York: Springer (2011). p. 43–62.

Google Scholar

55. Berg HC, Purcell EM. Physics of chemoreception. Biophys J (1977) 20(2):193–219. doi:10.1016/S0006-3495(77)85544-6

PubMed Abstract | CrossRef Full Text | Google Scholar

56. DeLisi C. The effect of cell size and receptor density on ligand-receptor reaction rate constants. Mol Immunol (1981) 18(6):507–11. doi:10.1016/0161-5890(81)90128-0

CrossRef Full Text | Google Scholar

57. Maity PC, Yang J, Klaesener K, Reth M. The nanoscale organization of the B lymphocyte membrane. Biochim Biophys Acta (2015) 1853(4):830–40. doi:10.1016/j.bbamcr.2014.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Yang J, Reth M. Oligomeric organization of the B-cell antigen receptor on resting cells. Nature (2010) 467(7314):465–9. doi:10.1038/nature09357

CrossRef Full Text | Google Scholar

59. Yang J, Reth M. The dissociation activation model of B cell antigen receptor triggering. FEBS Lett (2010) 584(24):4872–7. doi:10.1016/j.febslet.2010.09.045

CrossRef Full Text | Google Scholar

60. Perelson AS, DeLisi C. Receptor clustering on a cell surface. I. Theory of receptor cross-linking by ligands bearing two chemically identical functional groups. Math Biosci (1980) 48:71–110. doi:10.1016/0025-5564(80)90017-6

CrossRef Full Text | Google Scholar

61. Perelson AS, Weisbuch G. Immunology for physicists. Rev Mod Phys (1997) 69(4):1219–67. doi:10.1103/RevModPhys.69.1219

CrossRef Full Text | Google Scholar

62. Mech F, Wilson D, Lehnert T, Hube B, Thilo Figge M. Epithelial invasion outcompetes hypha development during Candida albicans infection as revealed by an image-based systems biology approach. Cytometry A (2014) 85(2):126–39. doi:10.1002/cyto.a.22418

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Medyukhina A, Timme S, Mokhtari Z, Figge MT. Image-based systems biology of infection. Cytometry A (2015) 87(6):462–70. doi:10.1002/cyto.a.22638

CrossRef Full Text | Google Scholar

64. Buhlmann D, Eberhardt HU, Medyukhina A, Prodinger WM, Figge MT, Zipfel PF, et al. FHR3 blocks C3d-mediated coactivation of human B cells. J Immunol (2016) 197(2):620–9. doi:10.4049/jimmunol.1600053

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Dustin ML. Signaling at neuro/immune synapses. J Clin Invest (2012) 122(4):1149–55. doi:10.1172/JCI58705

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Batista FD, Iber D, Neuberger MS. B cells acquire antigen from target cells after synapse formation. Nature (2001) 411(6836):489–94. doi:10.1038/35078099

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Weikl TR, Lipowsky R. Pattern formation during T-cell adhesion. Biophys J (2004) 87(6):3665–78. doi:10.1529/biophysj.104.045609

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Goodridge HS, Reyes CN, Becker CA, Katsumoto TR, Ma J, Wolf AJ, et al. Activation of the innate immune receptor Dectin-1 upon formation of a ‘phagocytic synapse’. Nature (2011) 472(7344):471–5. doi:10.1038/nature10071

CrossRef Full Text | Google Scholar

69. Ketchum C, Miller H, Song W, Upadhyaya A. Ligand mobility regulates B cell receptor clustering and signaling activation. Biophys J (2014) 106(1):26–36. doi:10.1016/j.bpj.2013.10.043

PubMed Abstract | CrossRef Full Text | Google Scholar