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.

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: (3)R+L⇄konmacrokoffmacro C. Here, konmacro is the reaction rate for binding, koffmacro is the dissociation rate and the corresponding association constant K_a is defined by their ratio: Ka=konmacro∕koffmacro.

The reaction equation (3) was then translated into the coupled system of ODE: (4)dRdt=−konmacroRL+koffmacroC, (5)dLdt=−konmacroRL+koffmacroC, (6)dCdt=+konmacroRL−koffmacroC. 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): (7)dCdt=αC2−βC+γ, where we defined the constants (8)α=konmacro, (9)β=koffmacro+konmacro [R(0)+L(0)], (10)γ=konmacro R(0)L(0). The ODE for C(t) can be solved by the separation of variables and yields the analytical solution: (11)C(t)=C− C+  1− eα(C+−C−)tC−− C+ eα(C+−C−)t with (12)C±=β2α±β24α2−γα . 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).

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

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 koffmicro=0.1 s−1, while the microscopic rate for RL binding was set to konmicro=106 s−1 (Figure 4A) and konmicro=107 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.

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 konmacro. 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: koffmacro=koffmicro. Arguments for this relation between macroscopic and microscopic dissociation rates are provided based on the analysis in Supplementary Material.

At the quantitative level, we observed differences in the binding kinetics depending on the receptor properties as well as on the microscopic binding rate konmicro. As could be expected, formation of RL complexes occurred slower for smaller konmicro=106 s−1 (Figure 4A) than for larger konmicro=107 s−1 (Figure 4B). Moreover, for a fixed value konmicro, 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 konmacro. The fitted curves are shown in Figure 4 and yielded for konmicro=106 s−1 (Figure 4A) the values konmacro≈1.9 µm3 s−1 for the ABM variants O-SOL and Y-SOL and konmacro≈0.6 µm3s−1 for the ABM variants O-MEM and Y-MEM. For konmicro=107 s−1 (Figure 4B), we obtained the values konmacro≈10.5 µm3s−1 for the ABM variants O-SOL and Y-SOL and konmacro≈1.7 µm3s−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 104 s−1≤konmicro≤106 s−1. Even though for konmicro>106 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.

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 koffmicro = koffmacro=0.1 s−1 and varying the microscopic binding rate in the range 104 s−1≤konmicro≤2.5×107 s−1. The corresponding macroscopic binding rate konmacro 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 konmacro(konmicro) 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.

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 konmacro(konmicro) in Figure 5 resemble Hill functions, (13)konmacrokonmicro=a konmicrob+konmicro, with parameters a and b that are specific for given receptor properties. Here, a denotes the upper limit for the macroscopic binding rate, konmacro(konmicro≫b)→a, and b is a constant that determines the slope of the Hill function, konmacro(konmicro≪b)→(a∕b)konmicro, while at intermediate value konmicro=b the Hill function attains half of its maximal value: konmacro(konmicro=b)=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 konmacro on konmicro 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 (14)konmacro(κ)=ks κks+κ, where ks=4π(rL+rR)(DL+DR) 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, κ=Vr konmicro, which is directly related to the microscopic binding rate konmicro and the reaction volume Vr=(4∕3)π(rL+rR)3 (35, 36). Combining equations (13) and (14) yields the following relationships: (15)a=ks, (16)b=ksVr. 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 (17)fs=aks, (18)fr=fs ksb Vr 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 fs 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 konmacro as a function of konmicro, 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 fr 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 fr Y−SOL of ABM variant Y-SOL was found to be similar to fr O−SOL with a relative decrease of only 4%, whereas this scaling factor for the ABM variants with membrane-anchored receptors, i.e., fr O−MEM and fr Y−MEM, was decreased by 74 and 61%, respectively. Furthermore, as shown in Figure 6, the scaling factors fs O−MEM and fs Y−MEM for membrane-anchored receptors were found to be decreased from fs 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.

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 konmacro as a function of konmicro. 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 konmacro on konmicro 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.

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 konmicro 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 konmicro=1×107 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 konmicro. 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 (NpR). The binding kinetics of ABM variants with NR=2×NpR monovalent receptors turned out to be even faster as the binding kinetics of bivalent ABM variants with NR=NpR (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.