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.


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 uniand 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)(15)(16)(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.

Microscopic Modeling of Receptor-Ligand Binding
Agent-based models (ABMs) are widely used in computational biology to simulate processes at the microscopic scale (9)(10)(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.

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)(21)(22)(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.

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   Table S1 in Supplementary Material. the specific distance ∆s = √ 2dD∆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 micro off and translates into the probability p micro off = k micro off ∆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 micro off = k macro off .
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 micro on = k micro on ∆t, where k micro on denotes the microscopic binding rate with unit s −1 . Note, that this rate is conceptually different from the macroscopic reaction rate k macro on with unit μm 3 s −1 , 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 micro on and k macro on by mapping the microscopic and macroscopic RL binding kinetics onto each other.

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 randomperiodic 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 from the smallest considered rate of binding (k micro on ) and dissociation (k micro off ) 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 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 5 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 9 . 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 −2 , i.e., reducing the B cell size by a factor 10 and the number of BCR to the order 10 3 . The ABM framework was implemented in the object-oriented programming language C++.

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: Here, k macro on is the reaction rate for binding, k macro off is the dissociation rate and the corresponding association constant K a is defined by their ratio: The reaction equation (3) was then translated into the coupled system of ODE: Assuming that initially no molecular complexes exist, where we defined the constants The ODE for C(t) can be solved by the separation of variables and yields the analytical solution: 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).

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 macro on and k micro on . 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 macro on (k micro on ) 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 macro on 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 macro on were used as data points to fit the optimal parameter values of the Hill equation k macro on (k micro on ) (see equation (13)) in order to map the microscopic and macroscopic binding kinetics.

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 −2 , 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.

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 micro off = 0.1 s −1 , while the microscopic rate for RL binding was set to k micro on = 10 6 s −1 ( Figure 4A) and k micro on = 10 7 s −1 (Figure 4B). Note that we provide the concentration of molecular complexes in units 1/μ m 3 to enable the comparison of the binding dynamics simulated by ODE and ABM variants with soluble and membraneanchored 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 4 molecules-yielded vanishingly small standard deviations.
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

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 micro on . As could be expected, formation of RL complexes occurred slower for smaller k micro on = 10 6 s −1 (Figure 4A) than for larger k micro on = 10 7 s −1 (Figure 4B). Moreover, for a fixed value k micro on , 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 membraneanchored 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 macro on . The fitted curves are shown in Figure 4 and yielded for k micro on = 10 6 s −1 ( Figure 4A) the values k macro on ≈ 1. Even though for k micro on > 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.

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  Table S7 in Supplementary Material for the four ABM variants.
The observed functional dependence of k macro on on k micro on is in agreement with theoretical considerations by Collins and Kimball on binding reactions of diffusing receptors and ligands in three spatial dimensions (31)(32)(33). They arrived at the expression where k s = 4π(r L + r R )(D L + D R ) denotes the diffusioncontrolled 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, κ = V r k micro on , which is directly related to the microscopic binding rate k micro on and the reaction volume V r = (4/3)π(r L + r R ) 3 (35,36). Combining equations (13) and (14) yields the following relationships: 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 r V 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 s k s with scaling factor f s . It then follows that f r and f s can be computed from the equations f r = 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. 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 macro on as a function of k micro on . This implies that the difference between ABM variants with soluble and membraneanchored 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 2 s −1 ) and membrane-anchored (D R = 0.05 μm 2 s −1 ) 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. 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.

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 micro on between 5 × 10 6 and 2.5 × 10 7 s −1 . 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 micro on = 1 × 10 7 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 micro on . 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.
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 (N R p ). The binding kinetics of ABM variants with N R = 2 × N R p monovalent receptors turned out to be even faster as the binding kinetics of bivalent ABM variants with N R = N R p (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.

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 micro off = k macro off = 0.1 s −1 the nonlinear relationship between the binding rates k macro on and k micro on . This resembles a Hill-type function (see Figure 5), which is in line with theoretical predictions by Collins and Kimball (31)(32)(33). Scanning  (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)(42)(43)(44)(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)(52)(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 membraneanchored 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)(58)(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 imagebased 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)(66)(67)(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.

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 micro on = 2.5 × 10 7 s −1 . The system is downscaled 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 −1 and the simulation time between two consecutive frames is 6.8 × 10 −8 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 micro on = 2.5 × 10 7 s −1 . The system is downscaled 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 −1 and the simulation time between two consecutive frames is 6.8 × 10 −8 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 micro on = 2.5 × 10 7 s −1 . 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 −1 and the simulation time between two consecutive frames is 6.8 × 10 −8 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 micro on = 2.5 × 10 7 s −1 . The system is downscaled 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 −1 and the simulation time between two consecutive frames is 6.8 × 10 −8 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 micro on = 2.5 × 10 7 s −1 . 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 −1 and the simulation time between two consecutive frames is 6.8 × 10 −8 s. A high-resolution video is available for download from https://asbdata.hki-jena.de/ LehnertFigge2017_FrontImmun/.