Mesoscale simulations of membrane-tethered reactions to parameterize cell-scale models of signaling

Volume vs. surface bi nding kinetics for a simple dimerization reaction. For bi molecular reactions, SpringSaLaD determines the microscopic probability of 2 binding sites forming a bond as they diffuse within a reaction radius that is slightly larger than the sum of their physical radii. The input to the algorithm is simply the macroscopic volumetric on rate constant (𝑘𝑘𝑜𝑜𝑜𝑜 (𝑣𝑣𝑜𝑜𝑣𝑣), units of µM-1s-1) and the diffusion coefficient of the individual spheres. Full details on the derivation of the reaction probability and a thorough validation of its accuracy can be found in the original paper describing SpringSaLaD (11). Importantly for the purposes of this work, the rate of binding for sites that happen to be tethered to a membrane are still treated as volumetric, because the spherical sites are located in the volume compartment even while they are constrained with links to the 2D membrane surface. .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint Figure 1 illustrates this for a simple dimerization reaction where the 2 yellow binding sites are tethered to the membrane anchor (gray sphere) by a 5nm link; in these simulations both the anchor and tether sites are given identical diffusion coefficients of 1µm2/s. We ran 100 SpringSaLaD simulations each with 40 dimerizing molecules (Fig. 1 only shows 2 for clarity). The mean trajectories for these 100 runs are then fit to a deterministic mass action membrane binding model in terms of surface densities using either COPASI or Virtual Cell (although an analytical solution can be also fit for simple dimerization) (Step 2 in Fig. 1). The output binding rate constant, 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚) in units of µm2molecule-1s-1, can then be used to parameterize larger deterministic or stochastic models with molecule numbers (>1000) or timescales (>10 s) that would be too large for even highly coarse-grained molecular simulators like SpringSaLaD. Also, we emphasize that binding rates are typically determined experimentally using in vitro volumetric measurements; that SpringSaLaD uses volumetric rate constants as its inputs makes the procedure in Fig. 1 especially appropriate and convenient. The results in Fig. 1 show that a 2nm diameter binding site with an on rate for dimerization of 0.47 µM -1s-1 and tethered to a membrane surface through a 5nm link, can be modeled as a 2D surface reaction with an on rate, 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚) , of 0.109 µm2molecule-1s-1. Using idealized models, we now explore how various structural and biophysical parameters control 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚)and when mass action rate constants may not be appropriate to describe membrane-associated binding kinetics. Figure 1. Workflow for finding membrane bimolecular binding rate constants, 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚), in terms of surface densities. Step 1 is to run 100 SpringSaLaD trajectories with the coarse grained model of the binding partners. Illustrated are a pair of simple monomers (top) consisting of 2nm diameter spheres (yellow) tethered to a 1nm membrane anchor sphere (gray) with a 5nm link; both the anchor and binding spheres are assigned the diffusion coefficient of 1 µm 2/s. At the bottom, the product dimer is depicted. Step 2 consists of fitting the average of 100 outputs from stochastic SpringSaLaD (SS) simulations to a deterministic (ODE) non-spatial model of surface-bound dimerization. = 0.47µM-1s-1 0 50 100 150 200 250 300 350 400 450 500 0 0.05 0.1 0.15 Surface Density of Dimer (molecule/µm2) Time (s) Fit = 0.109 µm2molecule-1s-1 Membrane Dimer (SS) Membrane Dimer (fit) Step 1: Run 100 dimerization Simulations in SpringSaLaD Step 2: Fit mean of the 100 SpringSaLaD trajectories to find in a deterministic dimerization model using surface densities .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint Effects of volumetric on-r ate, 2D diffusion rate, surface density and structural features on dimerization kinetics We s tart by analyzing the case of a 5nm stiff tether between a 1nm diameter binding sphere and the membrane. The diffusion coefficient for the binder sphere is set to 1µm2/s. Table 1 gives

for all combinations of two volumetric on rate constants, two surface densities and two membrane diffusion coefficients. All the dimerization rate laws are irreversible except for the last 2 rows, where the off-rate constant is indicated. The 2D on-rate constant derived by fitting the SpringSaLaD simulation, 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚), is in the 5th column; for consistency, all these fits are performed for kinetics at 80% completion. For comparison we also provide the 2D on-rate constant, 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ) = 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑣𝑣𝑜𝑜𝑣𝑣)/ℎ, where ℎ is the linker length plus the radius of the binder sphere. If ℎ is in units of µm, dividing by a unit conversion factor of 602.2 converts 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ) from units of µM-1s-1 to units of µm2molecules-1 s-1; 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ) is an equivalent 2D binding rate constant of freely diffusing monomers confined within a thin layer with a height ℎ adjacent to the membrane. The last column provides the ratio of the sum of the squared deviations to the sum of the squared SpringSaLaD mean values; this ratio provides a measure of the goodness of fit to the bimolecular mass action rate law, with anything less than ~10-3 representing a good fit. Examples of the fits are shown in Figure 2. We l ooked at two anchor diffusion coefficients corresponding to that of a large transmembrane protein domain (𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚= 0.01 µm2/s) and a lipid anchor (𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚=1 µm2/s). The first row of Table 1 considers a case where diffusion of the binder (𝐷𝐷𝑣𝑣𝑜𝑜𝑣𝑣=1 µm2/s) is 100 times that of the anchor site in the membrane. Because the linker is stiff, the anchor acts as a pivot and the binder rapidly moves within a hemispherical shell to effectively create a reaction region with a thickness of slightly larger than the 1-nm diameter of the binder sphere. Because the region of spatial overlap of the two shells where the binding may occur is restricted, it might seem surprising that 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚) is so closely approximated by 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ). However, as shown in the Supporting Information, the reaction probability is enhanced because of the effectively higher density of binding sites .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint within this restricted region, which compensates for the smaller spatial overlap. Thus, the calculations in the Supporting Information both explain and validate the SpringSaLaD results. The second row of Table 1 corresponds to the case where the anchor and the binder have the same fast diffusion. In this scenario, that 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚)≅ 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ) is intuitive, because the effect of the tether in this case essentially reduces to confining the binders within the layer adjacent to the membrane. The solutions in the Supporting Information assume that the molecular distributions are spatially uniform at any time, which pertains to the cases of the first two rows of Table 1 where 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ)<𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚. I t has been shown that the rates of membrane reactions may be susceptible to deviations from a simple mass action rate law (13-17), which manifests themselves as significant changes of the apparent rate constant with initial surface density. These changes need to be assessed before using the rate constants we obtain by the procedure of Fig. 1 in large scale cell-level models. We probed for this by decreasing the initial density by a factor of 100. The combinations of 𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚and 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑣𝑣𝑜𝑜 𝑣𝑣) in the third and fourth rows of T able 1 resulted in relatively small changes in 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚), indicating that the mass action rate law applies to these cases. We further tested this by increasing the 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑣𝑣𝑜𝑜 𝑣𝑣) by a factor of 100 in the lower half of Table 1. Clearly, for the case of 𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚= 0.01 µm2/s, there is a strong dependence on surface density. Furthermore, this value of 𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚 is below the “well-mixed” limit given by 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ). Thus, the combination of 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑣𝑣𝑜𝑜 𝑣𝑣)= 1.0 µM-1s-1 and 𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚= 0.01 µm2/s present cases where mass action rate laws would not apply. In general, it would not be appropriate to use a mass action rate law, if the apparent 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚) is significantly greater than 𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚 . Of course, the appropriateness of the mass action rate law can also be judged by the goodness of fit when determining 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚) from the SpringSaLaD simulation (last column of Table 1); as demonstrated in Fig. 2, the case of 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑣𝑣𝑜𝑜 𝑣𝑣)= 1.0 µM-1s-1 and 𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚= 0.01 µm2/s is not fit very well to a mass-action dimerization rate. The simulation results are initially faster and then ultimately slower than the best fit that assumes mass action kinetics. This is because the monomers whose binding sites are initially close to each other (effectively within the “reach” of the tether) will react, but leave behind depletion zones where monomers are too far away from potential binding partners (15, 16). These isolated monomers can be discerned toward the end of Movie 1, which presents an example trajectory for this case. Figure 2. Examples of SpringSaLaD simulation data and their fit to surface-confined mass action dimerization kinetics. The initial surface density of monomers in each case is 2500 molecule/µm2. The diffusion coefficient of the anchor and the volumetric binding rate constant of the binding sites are indicated above each graph, corresponding respectively to rows 1, 6 and 5 of Table 1. Each of the SpringSaLaD simulations were run to 80% completion. .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint Consistent with the idea of depletion zones, the last 2 rows of Table 1 show that when reversibility is introduced, the fitted on-rate constant increases (compare, respectively, to rows 5 and 7 of Table 1). This is because when dimerization is reversible, free monomers can reappear to fill in depletion zones, thereby countering the slow diffusion. While the fits to mass action kinetics are still poor, especially for the low density case, the respective mean steady-state densities of dimers in the SpringSaLaD simulations, 570 and 5.7 molecules/µm2, are consistent with the thermodynamic law of mass action (18). According to this law, the ratio of the squared monomer volumetric concentration and the dimer volumetric concentration is determined at steady state by the dissociation constant, 𝐾𝐾𝑑𝑑= 𝑘𝑘𝑜𝑜𝑜𝑜𝑜𝑜/𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ). Tabl e 2 provides results for three computational experiments in which structural features of the SpringSaLaD molecules are varied. These are all for dimerization reactions where the maximum distance between the membrane and the binding site is 4 times longer than in Table 1: ℎ = 0.0205 µm (linker length of 20nm and binding site radius of 0.5nm). The 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ) is therefore a factor of 4 slower than that in Table 1 for the same 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑣𝑣𝑜𝑜𝑣𝑣) of 1µM-1s-1. The SpringSaLaD simulations were carried out with a slow membrane-anchor diffusion coefficient, 𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚= 0.01µm2/s, to model a membrane receptor. The first row can be directly compared to the fifth row of Table 1, where the only difference is the length of the linker. As would be expected from the increase in h, the membrane on rate constant is decreased; however, importantly, this constant is now closer to the diffusion rate and therefore deviations from mass action are significantly reduced. The second row in Table 2 shows results for a structure where additional spherical sites are introduced between the membrane anchor and the binding site to model the space occupied by a cytosolic protein sequence; this steric effect results in a small decrease in 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚). This decrease is reversed when flexibility is introduced by allowing the spherical sites to be pivot points in the third row of the Table; this is how disordered domains may be modeled in SpringSaLaD. Overall, for these idealized structures and simple dimerization, the fitted on-rate constants in Table 2 are relatively close to 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ). Application of the method to interaction of receptor-bound SOS with Ras. Till now, we have employed idealized molecular structures to validate our method and to learn some biophysical principles that control the on-rate constants of binding sites tethered to membranes. We now illustrate the application of this approach to a biologically relevant .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint example, namely the interaction of the G-protein exchange factor (GEF) SOS with the lipid anchored small G-protein Ras (19). SOS has 2 binding sites for Ras: an allosteric site and a catalytic site. When a Ras molecule binds to the allosteric site it increases the GEF activity of the catalytic site (20). Additionally, before SOS binds to Ras it is first recruited to an active receptor tyrosine kinase (RTK) through an adapter protein; the adaptor binds to a proline rich motif (PRM) on SOS via a SH3 domain and to a phosphorylated tyrosine via a SH2 domain. One such RTK is the Epidermal Growth Factor Receptor (EGFR) and one such adaptor protein is Grb2 (21). Once SOS is bound to Grb2, it becomes membrane tethered and its interaction with Ras is facilitated (20). However the complex mechanistic details are still emerging (22). Figure 3. Membrane binding of Ras to the catalytic site of Receptor-bound SOS. Top: direct binding. Bottom: after pre-binding of Ras to the allosteric site. The molecular structures are approximated by using the SpringSaLaD 3D editing utility based on atomic structures derived from AlphaFold2. The top left structure is an EGFR cytoplasmic domain anchored to the membrane (red kinase domain, followed by yellow disordered tail capped by a phosphotyrosine in green); the latter is linked to a cyan SH2 domain in Grb2; one of its magenta SH3 domains is linked to an olive PRM on the end of the disordered region of SOS; the violet SOS binding site for the catalytic domain of Ras is indicated with an asterisk (*). The bottom left structure is identical, except that the pink allosteric site on SOS is bound to Ras. The Ras structures are shown in the center with the yellow binding sites indicated by an asterisk. The input rate constants for the SpringSaLaD simulations are shown at the top, corresponding to the volumetric on rate for Ras binding to the catalytic site of SOS. The EGFR anchor diffusion coefficient is 0.01µm 2/s. All other site diffusion coefficients are 1.0µm2/s. For each condition, 20 EGFR-Grb2-SOS molecules react with 200 Ras molecules on a 250nmX250nm membrane surface to generate 100 SpringSaLaD trajectories. Their means were fitted to a deterministic 2D rate law to derive 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚), fixing 𝑘𝑘𝑜𝑜𝑜𝑜𝑜𝑜 at 4.0 s-1; results for the 2 conditions are shown on the right. .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint We asked the limited question of how binding of Ras with the receptor-associated SOS catalytic site might depend on whether SOS is prebound to Ras at the allosteric site. Ras is a lipid- anchored protein, so we reasoned that binding of Ras to the allosteric site of SOS would bring the SOS catalytic site closer to the membrane to enhance binding to a second Ras and subsequent exchange of GDP for GTP . Just how big an effect this is, may be estimated by the procedure developed above, with the results shown in Fig. 3. We dev eloped molecular models with the aid of the mol2sphere (23) utility within SpringSaLaD and were guided by AlphaFold 2 atomic structure predictions (24, 25); all the site diameters and linker lengths are available in the SpringSaLaD input file included in the supporting information; snapshots of the structure are shown in Fig. 3. The top of Fig. 3 displays results for binding of SOS-Grb2-EGFR to Ras at the SOS catalytic site; the bottom shows results for the same reaction, except SOS-Grb2-EGFR had been first bound to a Ras molecule at the SOS allosteric site. The on and off rates shown at the top of Fig. 3 (19) are applied to both of the reactions. In these models, the EGFR membrane anchor site is assigned a diffusion coefficient of 0.01 µm 2/s to represent a large transmembrane protein, while the Ras membrane anchor is assigned a diffusion coefficient of 1.0 µm 2/s to represent a lipid anchor; all the sites that are dangling in the cytosol volume are given 𝐷𝐷𝑣𝑣𝑜𝑜 𝑣𝑣 of 1.0 µm2/s, but since the binding reaction is not close to diffusion-limited, the precise values are not critical. Consistent with these being reaction-limited on-rates, the SpringSaLaD simulation outputs (averages of 40 runs) are well fitted to reversible mass action kinetic law, as shown in the plots on the right of Fig. 3 (relative squared deviations are, respectively, 4.1 X 10-4 and 2.0 X 10-4). Importantly, the 2D on-rate constants (𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚)) derived from these fits are, respectively, 1.1 X 10-3 µm2molecules-1s-1 and 4.0 X 10-3 µm2molecules-1s-1. Likewise, the affinity of the catalytic is site is increased by allosteric site pre- association: 𝐾𝐾𝑑𝑑 = 3600 molecules/µm2 for the top of Fig. 3 and 1000 molecules/µm2 for the bottom pre-association case. Thus, SOS allosteric site association with Ras is estimated to accelerate its catalytic site binding and affinity by a factor of ~4 – even when SOS is already confined to the membrane through Grb2-mediated association with EGFR. D iscussion The k inetics of reactions at membranes have long fascinated biophysicists (1, 3-5, 8, 10, 14, 17). These studies have produced theoretical insights to illuminate how surface-associated reactions have distinct properties compared with reactions occurring in 3D solution. Which of these special properties are most pertinent to any given membrane-bound molecular interaction is difficult to ascertain a priori. Furthermore, experiments to measure bimolecular kinetics on membrane surfaces are difficult, so often only on-rate constants measured in 3D are accessible. Fundamentally, however, the kinetics of key membrane-associated reactions depend on the 2D surface densities and 2D rate constants, not on the bulk cellular concentrations and 3D rate constants. Indeed, because surface to volume ratios of different cell types vary tremendously, volumetric rate constants cannot be readily used to model and simulate cell signaling systems. To address these theoretical and practical problems, we describe a procedure (Fig. 1) using experimentally accessible volumetric on-rate constants, 𝑘𝑘 𝑜𝑜𝑜𝑜 (𝑣𝑣𝑜𝑜 𝑣𝑣) with the SpringSaLaD simulation software to estimate the 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚), the 2-dimensional rate constant for a membrane-confined bimolecular reaction. T o validate the method, we applied it to the dimerization of a single binding site tethered to a surface through a 5nm stiff linker, where the membrane anchor acts as a pivot (Table 1). For the situation where the reaction is rate limiting, this system can be solved analytically (see .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint Supporting Information); gratifyingly, 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚) determined by our method is well reproduced by the analytical solution. Interestingly, for these cases, 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚) is well approximated by 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ) = 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑣𝑣𝑜𝑜 𝑣𝑣)/ℎ/602.2 (µm2molecules-1s-1), where ℎ is the distance of the binding site from the membrane anchor. This parameter has also been referred to as the “confinement length” (8), defining a thin volume above the membrane that concentrates the binding sites and directly producing the relationship between 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ) and 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑣𝑣 𝑜𝑜𝑣𝑣). W hile the mass action kinetics are generally applicable for both encounter-limited and reaction- limited kinetics in 3D solution (but see (15)), it has long been appreciated that the situation may be more complex for 2D kinetics (3, 5, 14, 17, 26). This is demonstrated by the results in Table 1 for situations where the diffusion coefficient of the anchor is slow, but the volumetric on-rate constant is fast. For these cases, different estimates of 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚) are obtained at different initial surface density – clearly incompatible with the mass action kinetics. Indeed, the third panel of Figure 2 shows that the SpringSaLaD kinetic data is not well fitted by a 2D mass action rate law. A video of one of these trajectories (Movie 1) nicely illustrates how the initial rate is fast, while the binding sites are within “reach” (10), but falls off as binding sites are left orphaned outside the reach of the remaining slowly diffusing monomers. The results allow us to generalize that mass action applies as long as 𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚 is close to or greater than 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ). Furthermore, a mass action rate law applies well to the initial rate, before depletion zones develop. Thus, our analysis provides an approach to determine whether 2D binding might be well approximated by the mass action kinetics, and if so, to estimate the 2D mass action on-rate constant. To explore how other molecular structural features might affect dimerization of the monomers tethered to the membrane, we looked at 3 additional idealized systems in Table 2. In all these, h was 20.5nm (as opposed to 5.5nm in Table 1). As expected, the longer confinement length decreased the estimated 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚) by a factor of about 4. The insertion of steric sites between the anchor and the binding site or allowing for flexibility of the linker region have minor effects on 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚), which is relatively well approximated by 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ) . Importantly, the lower value of 𝑘𝑘𝑜𝑜𝑜𝑜 (ℎ) for this system was closer to 𝐷𝐷𝑚𝑚𝑚𝑚𝑚𝑚, resulting in a better fit by the mass action kinetics than for the similar case with the 5nm linker (5th row of Table 1). For membrane binding of real biological molecules, it may be difficult to estimate the average location of a binding sites relative to the membrane surface (i.e. h); also, the 2 binding sites may be parts of very different structures with different distances from the membrane surface. In situations like this, our approach has the potential to provide good estimates of rate constants that can be applied to larger cell signaling systems. Indeed, there may be direct insights that can be realized just by considering the structural details of the interacting membrane molecules. We have illustrated this in relation to adaptor-mediated protein kinase receptor signaling mechanisms, specifically for the interaction of the GEF SOS with its effector Ras (Fig. 3). It has been shown that direct catalysis by the SOS catalytic domain of Ras conversion from the GDP to the GTP states is relatively slow. However, prebinding binding of Ras to SOS at a site that is not catalytic (termed the “allosteric” site on SOS) significantly accelerates the catalytic activity, where the catalysis becomes processive (19, 20, 22). The results in Fig. 3 suggest that at least part of this acceleration may be due to the close proximity of the SOS catalytic site to the membrane once it is bound to Ras at its allosteric site. Even though SOS is already localized to the membrane by initially binding to EGFR via Grb2 in our computational experiment, pre- binding of the SOS allosteric site to Ras brings it to still closer proximity to the membrane. Of course, there could be additional effects such as a direct allosteric enhancement through a conformational change or release of self-inhibition (19, 20, 22), but here we focus on the .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint significance of constricting the binding zone through membrane tethers of varying length and flexibility. Our approach toward deriving membrane on-rates will aid in the parametrization of ODE and PDE models that could help elucidate the full kinetic and mechanistic details of processes such as Ras activation and downstream signaling.

All simulations were performed with SpringSaLaD v. 2.3.4 (https://vcell.org/ssalad). To build coarse grained molecular structures used in SpringSaLaD for Fig. 3, we used AlphaFold2 (24, 25) to generate PDB file estimates of protein structures for EGFR, Grb2, SOS, and Ras via input of entire amino acid sequences. These PDB files are converted to coarse-grained molecular models via the mol2sphere (23) utility embedded in SpringSaLaD. In some cases, we manually edited the structures to capture their essential features from measurements on the PDB structures, as visualized in PyMol (Schrödinger, Inc.). In particular, for SOS we subdivide the CDC25 and REM domains into multiple smaller spherical sites with only one site capable of participating in a binding reaction. This process maintains the structural characteristics of these domains, while ensuring that the binding radius of the domain is not artificially inflated. We also consider molecule flexibility when making user modifications to molecule structure. When PDB files are imported to SpringSaLaD via mol2sphere, the default is for each site to be linked to no more than 2 immediately adjacent sites. Flexibility can be decreased by introducing more stiff linkers to more adjacent spherical sites. We employ this method to ensure that our multi-site representation of CDC25 and REM domains diffuse as a fixed group of spherical sites instead of individual, highly flexible domains. Modeling disordered regions, such as the PRM region of SOS, can be challenging due to low confidence in the AlphaFold2-generated geometry of these regions. To model disordered domains, we use PyMOL to measure the length of entire straight chain amino acid sequences, then model this sequence in SpringSaLaD using 1.0 nm diameter sites connected by 3.1 nm linkers. Binding reactions in all simulations have rates input in terms of µM -1s-1; successful binding results in 1nm (Tables 1 and 2) or 0.5 nm (Fig. 3) distances between the surfaces of the spherical sites. A set of 100 trajectories are simulated in parallel using the Center for Cell Analysis and Modeling High Performance Compute Cluster (https://health.uconn.edu/high-per formance-computing/resources/). SpringSaLaD input files are in Supporting Information and provide all the geometric details for the molecules in each computational experiment. The m ean of 100 SpringSaLaD simulations for each of the simulations in Tables 1 and 2 and Fig. 3 were fit to a deterministic 2D binding rate law to obtain 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚). For irrreversible dimerization (Table 1, first 8 rows and Table 2), a fit to an analytical expression (Eq. 1) used the Excel solver; for reversible dimerization (Table 1, rows 9 and 10) and for the fits in Fig. 3, we used the COPASI (27) parameter estimation tool within Virtual Cell (VCell) (28, 29). The latter can be accessed in the VCell published BioModel “Peterson Figure 3: Ras- SOS_Binding_fit_to_SpringSaLaD”. All these results with some further analysis can also be found in the spreadsheets included in the Supporting Information. Acknowledgments This work was supported by NIH grants R24 GM137787 and R01 GM132859. We are pleased to acknowledge the advice of Aniruddha Chattaraj with some of the data analysis. .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint 1. J. M. Haugh, A Unified Model for Signal Transduction Reactions in Cellular Membranes. Biophys. J. 82, 591-604 (2002). 2. B. M. Slepchenko, J. C. Schaff, I. Macara, L. M. Loew, Quantitative cell biology with the Virtual Cell. Trends Cell Biol 13, 570-576 (2003). 3. G. Adam, M. Delbrück, "Reduction of Dimensionality in Biological Diffusion" in Structural Chemistry and Molecular Biology, A. Rich, N. Davidson, Eds. (W. H. Freeman and Co., San Francisco, 1968), pp. 198-215. 4. M. A. McCloskey, M. Poo, Rates of Membrane-associated Reactions: Reduction of Dimentionality Revisited. Journal of Cell Biology 102, 88-96 (1986). 5. D. Axelrod, M. D. Wang, Reduction-of-dimensionality kinetics at reaction-limited cell surface receptors. Biophys J 66, 588-600 (1994). 6. B. Windisch, D. Bray, T. Duke, Balls and Chains—A Mesoscopic Approach to Tethered Protein Domains. Biophysical Journal 91, 2383-2392 (2006). 7. G. I. Bell, M. Dembo, P . Bongrand, Cell adhesion. Competition between nonspecific repulsion and specific bonding. Biophys J 45, 1051-1064 (1984). 8. Y . Wu, J. Vendome, L. Shapiro, A. Ben-Shaul, B. Honig, Transforming binding affinities from three dimensions to two with application to cadherin clustering. Nature 475, 510- 513 (2011). 9. Z.-R. Xie, J. Chen, Y . Wu, Linking 3D and 2D binding kinetics of membrane proteins by multiscale simulations. Protein Science 23, 1789-1799 (2014). 10. Y . Zhang et al., The Influence of Molecular Reach and Diffusivity on the Efficacy of Membrane-Confined Reactions. Biophys J 117, 1189-1201 (2019). 11. P . J. Michalski, L. M. Loew, SpringSaLaD: A Spatial, Particle-Based Biochemical Simulation Platform with Excluded Volume. Biophys J 110, 523-529 (2016). 12. A. Chattaraj, M. Youngstrom, L. M. Loew, The interplay of structural and cellular biophysics controls clustering of multivalent molecules. bioRxiv (accepted in Biophys. J.) 10.1101/373084, 373084 (2018). 13. R. Kopelman, Rate processes on fractals: Theory, simulations, and experiments. Journal of Statistical Physics 42, 185-200 (1986). 14. H. Berry, Monte Carlo Simulations of Enzyme Reactions in Two Dimensions: Fractal Kinetics and Spatial Segregation. Biophysical Journal 83, 1891-1901 (2002). 15. A. A. Ovchinnikov, Y . B. Zeldovich, Role of density fluctuations in bimolecular reaction kinetics. Chemical Physics 28, 215-218 (1977). 16. D. Toussaint, F. Wilczek, Particle–antiparticle annihilation in diffusive motion. The Journal of Chemical Physics 78, 2642-2647 (1983). 17. D. C. Torney, H. M. McConnell, G. R. Porter, Diffusion-limited reaction rate theory for two-dimensional systems. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 387, 147-170 (1983). 18. A. B. Koudriavtsev, R. F. Jameson, W. Linert, SpringerLink, The Law of Mass Action, Engineering Online Library (Springer Berlin Heidelberg : Imprint: Springer, Berlin, Heidelberg, ed. 1st 2001., 2001). 19. L. Iversen et al., Ras activation by SOS: Allosteric regulation by altered fluctuation dynamics. Science 345, 50-54 (2014). 20. P . Bandaru, Y . Kondo, J. Kuriyan, The Interdependent Activation of Son-of-Sevenless and Ras. Cold Spring Harb Perspect Med 9 (2019). 21. R. N. Jorissen et al., Epidermal growth factor receptor: mechanisms of activation and signalling. Experimental Cell Research 284, 31-53 (2003). 22. H. Ren, A. A. Lee, L. J. N. Lew, J. B. DeGrandchamp, J. T. Groves, Positive feedback in Ras activation by full-length SOS arises from autoinhibition release mechanism. Biophysical Journal 123, 3295-3303 (2024). .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint 23. J. Masison, P . J. Michalski, L. M. Loew, A. D. Schuyler, mol2sphere: spherical decomposition of multi-domain molecules for visualization and coarse grained spatial modeling. Bioinformatics 34, 3948-3950 (2018). 24. J. Abramson et al., Accurate structure prediction of biomolecular interactions with AlphaFold 3. Nature 630, 493-500 (2024). 25. J. Jumper et al., Highly accurate protein structure prediction with AlphaFold. Nature 596, 583-589 (2021). 26. G. I. Bell, Models for the Specific Adhesion of Cells to Cells. Science 200, 618-627 (1978). 27. S. Hoops et al., COPASI: a COmplex PAthway SImulator. Bioinformatics 22, 3067-3074 (2006). 28. J. Schaff, C. C. Fink, B. Slepchenko, J. H. Carson, L. M. Loew, A general computational framework for modeling cellular structure and function. Biophys J 73, 1135-1146 (1997). 29. A. E. Cowan, Moraru, II, J. C. Schaff, B. M. Slepchenko, L. M. Loew, Spatial modeling of cell signaling networks. Methods Cell Biol 110, 195-221 (2012). Figure Captions Figure 1. Workflow for finding membrane bimolecular binding rate constants (𝑘𝑘 𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚)) in terms of surface densities. Step 1 is to run 100 SpringSaLaD trajectories with the coarse grained model of the binding partners. Illustrated are a pair of simple monomers (top) consisting of 2nm diameter spheres (yellow) tethered to a 1nm membrane anchor sphere (gray) with a 5nm link; at the bottom, the product dimer is depicted. Step 2 consists of fitting the average of 100 outputs from stochastic SpringSaLaD simulations to a deterministic (ODE) non-spatial model of surface-bound dimerization. Figure 2. Examples of SpringSaLaD simulation data and their fit to surface-confined mass action dimerization kinetics. The initial surface density of monomers in each case is 2500 molecule/µm 2. The diffusion coefficient of the anchor and the volumetric binding rate constant of the binding sites are indicated above each graph, corresponding respectively to rows 1, 6 and 5 of Table 1. Figure 3. Membrane binding of Ras to the catalytic site of Receptor-bound SOS. Top: direct binding. Bottom: after pre-binding of Ras to the allosteric site. The molecular structures are approximated by using the SpringSaLaD 3D editing utility based on atomic structures derived from AlphaFold2. The top left structure is an EGFR cytoplasmic domain anchored to the membrane (red kinase domain, followed by yellow disordered tail capped by a phosphotyrosine in green); the latter is linked to a cyan SH2 domain in Grb2; one of its magenta SH3 domains is linked to an olive PRM on the end of the disordered region of SOS; the violet SOS binding site for the catalytic domain of Ras is indicated with an asterisk (*). The bottom left structure is identical, except that the pink allosteric site on SOS is bound to Ras. The Ras structures are shown in the center with the yellow binding sites indicated by an asterisk. The input rate constants for the SpringSaLaD simulations are shown at the top, corresponding to the volumetric on rate for Ras binding to the catalytic site of SOS. The EGFR anchor diffusion coefficient is 0.01µm 2/s. All other site diffusion coefficients are 1.0µm2/s. For each condition, 20 EGFR-Grb2-SOS molecules react with 200 Ras molecules on a 250nmX250nm membrane .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint surface to generate 100 SpringSaLaD trajectories. Their means were fit to a deterministic 2D rate law to derive 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚), fixing koff at 4.0 s-1; results for the 2 conditions are shown on the right. List of Supporting Data and Information Derivations for limiting cases related to Table 1. Filename: Supporting Text - Derivations for limiting cases.docx SpringSaLaD input files for the models in Tables 1 and 2 and Figure 3. Filenames: Table 1 and 2 SIMS.zip; EGFR_Grb_SOS binding to Ras.txt; EGFR_Grb_SOS_prebound at allo binding to Ras.txt Spreadsheets related to Tables 1 and 2 and Figure 3 containing the SpringSaLaD simulation

and the fits to determine 𝑘𝑘𝑜𝑜𝑜𝑜 (𝑚𝑚𝑚𝑚𝑚𝑚). Filenames: Table 1 5nm_stiff_fits.xlsx; Table2 20nmT_fits.xlsx; Data for Figure 3.xlsx Movie of an example trajectory for the model of Table 1 Row 5 and also Fig. 2, 3rd panel. Filename: Movie kon=1 Danchor=membrane protein.mp4 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted December 5, 2024. ; https://doi.org/10.1101/2024.12.04.626844doi: bioRxiv preprint