Computational framework for single-cell spatiotemporal dynamics of optogenetic membrane recruitment

Summary We describe a modular computational framework for analyzing cell-wide spatiotemporal signaling dynamics in single-cell microscopy experiments that accounts for the experiment-specific geometric and diffractive complexities that arise from heterogeneous cell morphologies and optical instrumentation. Inputs are unique cell geometries and protein concentrations derived from confocal stacks and spatiotemporally varying environmental stimuli. After simulating the system with a model of choice, the output is convolved with the microscope point-spread function for direct comparison with the observable image. We experimentally validate this approach in single cells with BcLOV4, an optogenetic membrane recruitment system for versatile control over cell signaling, using a three-dimensional non-linear finite element model with all parameters experimentally derived. The simulations recapitulate observed subcellular and cell-to-cell variability in BcLOV4 signaling, allowing for inter-experimental differences of cellular and instrumentation origins to be elucidated and resolved for improved interpretive robustness. This single-cell approach will enhance optogenetics and spatiotemporally resolved signaling studies.

ivan.kuznetsov@pennmedicine.upenn. edu (I.A.K.), bchow@seas.upenn.edu (B.Y.C.) In brief Kuznetsov et al. describe a computational framework for analyzing cell signaling in single-cell microscopy and use it to study BcLOV4 optogenetic membrane recruitment. By accounting for inter-experimental differences of cell and instrumentation origin, the approach allows for robustness to inter-experimental variability. This method will enhance spatiotemporally resolved signaling studies.

INTRODUCTION
Quantitative understanding of subcellular signaling dynamics is key to cell physiology and predicting cell behaviors. To this end, a growing number of increasingly complex models and model design tools have been reported, e.g., Virtual Cell (VCell), MCell/CellBlender, and COMSOL packages, among others (Stiles et al., 1996;Schutter, 2000;Kerr et al., 2008;Moraru et al., 2008;Menshykau, 2013;Hallett et al., 2016). While single-cell analyses can reveal the physiological underpinnings of heterogeneous cell behaviors, cell models often do not accurately reflect such heterogeneity. They typically are fit to popula-tion-average phenotypes and apply simplifying assumptions regarding cell geometry and protein expression (Cheong et al., 2010;Handly et al., 2016). Thus, validated computational frameworks to predict spatiotemporally resolved subcellular signaling dynamics of unique single cells can advance quantitative biology. Accounting for single-cell experimental variability from extrinsic or non-cellular factors, such as instrumentation-dependent diffractive phenomena in an optogenetics experiment, would further improve comparisons between model and experiment.
To these ends, we created a modular in silico framework for spatiotemporally resolved study of single-cell signaling in non-equilibrium dynamics across a unique three-dimensional (3D) cell and mirrors optical hardware-dependent experimental conditions used to generate single-cell data. Thus, it will enhance our biophysical understanding of signaling processes and guide the onward design of spatiotemporally complex signaling outputs.
We apply this framework to study dynamic membrane recruitment of a cytosol-sequestered protein, a ubiquitous approach in optogenetics and related chemical techniques to control signaling (Inoue et al., 2005;Repina et al., 2017;Stanton et al., 2018;Hannanta-Anan et al., 2019). It can be mediated by an inducible heterodimerizing protein-protein interaction or by inducible protein-lipid interaction with the plasma membrane itself in a single-component system (He et al., 2017;Glantz et al., 2018;Li et al., 2021). Here, we study the single-component BcLOV4 that we have previously used for versatile blue light control over GTPase signaling (Hannanta-Anan et al., 2019;Berlew et al., 2020Berlew et al., , 2021Berlew et al., , 2022. We employ a non-linear 3D finite element model (FEM). This is well suited to the complex geometries of mammalian cells and organelles, which are inherently highly irregular in morphology, size, protein expression, and protein distribution. It is unlike existing mathematical models of optogenetic membrane recruitment that have been idealized (symmetric twodimensional [2D]), have been solved under steady-state conditions and/or have assumed linearity with respect to concentration, despite the fact that the 3D recruitment process is dynamic and increasingly non-linear with membrane binding site occupancy (Valon et al., 2015;Niu et al., 2016;de Beco et al., 2018;Natwick and Collins, 2021). In brief, model biophysical parameters were experimentally derived. For each cell, the mesh was autogenerated from cell images and then seeded with the initial protein concentration. Use of actual cell geometry and protein concentration ensured accurate prediction of the spatiotemporal signaling outputs in response to arbitrary light inputs. Furthermore, inclusion of instrumentdependent diffractive effects, by way of PSF correction of the simulated output, permitted direct correlation with its corresponding single-cell image data.
We use the framework to explore how cell-specific morphology, optogenetic protein-specific biophysical properties, and the hardware-specific optical properties of the stimulatory input all shape the resultant observed signaling output. These analyses successfully identified the origins of detection artifacts that confound data interpretation. To demonstrate generality, we also show how existing software tools can be integrated into the modular framework by implementing it in VCell to similar effect. This work establishes a foundation for further spatiotemporally accurate models of peripheral membrane protein signaling and shows how single-cell models can reflect and combat experimental heterogeneity for improved interpretive robustness.

RESULTS
BcLOV4 is a well-posed and simple test bed for validating the framework. Its optical stimulation can be modulated with subcellular precision, and the diffractive properties of the optical hard-ware are measurable when testing how well the framework accounts for cell-extrinsic and experiment-specific factors. As BcLOV4 directly binds the plasma membrane in a light-dependent manner, its minimal single-component nature simplifies modeling by limiting the set of equations and the terms/ constants to be measured; it also streamlines experiments without the need to express and quantify multiple proteins in live cells. Henceforth, ''BcLOV4'' refers to the fluorescently labeled BcLOV4-mCherry, unless stated otherwise; also, ''lit state'' or ''lit'' refers to its photoactive state, and ''dark state '' or ''dark'' refers to dark-adapted protein. There are several components of the framework implementation described herein: (1) the automated mesh generation from image volumes of corresponding cells; (2) the set of non-linear partial differential equations (PDEs) and associated custom solver; (3) the experimental measurement of biophysical constants in the model; and (4) the convolution of the microscope PSF to correlate the theoretical output to the experimental observable.

Automated 3D mesh generation
The 3D mesh was necessary for biophysical accuracy. Since the stimulation volume has a non-negligible axial extent, the protein diffuses in three dimensions in and out of the imaging plane, to and from the membrane, and within the bulk cytosol ( Figure 1B). Diffractive effects lead to inevitable contributions of cytosolic fluorescence to membrane fluorescence signal and vice versa ( Figure 1B), which can confound data interpretation. This led us to generate a unique tetrahedral mesh for each experimental cell and construct a 3D FEM. Mesh generation was almost fully automated from fluorescence micrographs of HEK cells expressing BcLOV4. The shape of the plasma and nuclear membranes were reconstructed from confocal microscopy image stacks by interpolation in the z axis ( Figure 1C), or were extrapolated by hemi-ellipsoid projections (i.e., adherent cells on flat-bottom substrates) when only one focal plane was available ( Figure 1D). Mesh nodes were initialized with the corresponding protein concentration in the starting dark-adapted state.
The v M state includes the period that BcLOV4 stays associated with the membrane ($60 s) after thermal reversion of the flavin photoadduct (Glantz et al., 2018;Glantz et al., 2019). The set of PDEs defines a non-linear reaction-diffusion model dependent on the photoconversion, lipid binding kinetics, and diffusion rates of BcLOV4 (Equation 1). The nonlinearity reflects the decrease in available membrane binding sites with protein translocation. The PDEs specifically belong to a category of models termed bulk-surface models, due to the inherent protein coupling between the 3D cytoplasm (the bulk) and 2D membrane (the surface) (Equation 2). This maturing class of models is of paramount importance to studying cell signaling that occurs at the inner leaflet Cell Reports Methods 2, 100245, July 18, 2022 3 Article ll OPEN ACCESS (Madzvamuse and Chung, 2016;Cusseddu et al., 2019;Paquin-Lefebvre et al., 2020).
Bulk-surface FEM on stationary volumes has been described (Madzvamuse and Chung, 2016;Cusseddu et al., 2019;Paquin-Lefebvre et al., 2020). Here, we adapt the model equations to describe the given bulk-surface reaction-diffusion model for a unique 3D volume. Consider the cytoplasm as a volume U in R 3 enclosed by the plasma membrane, a closed hypersurface dU. We define u C ðx; tÞ; v C ðx; tÞ : U/R 3 as the cytosolic lit and dark-state protein, respectively, and similarly define u M ðx; tÞ; v M ðx; tÞ : dU/R 2 as the membrane-bound lit and dark-state protein, respectively; here, the units are molecules/ mm 3 for u C ; v C , and molecules/mm 2 for u M ; v M . We then pose the system of coupled PDEs: with the bulk-surface coupling enforced using Robin-type boundary conditions (Cusseddu et al., 2019): where D dU is the Laplace-Beltrami operator, n is the outwardfacing unit normal vector to U, D cyto (mm 2 /s) is the cytoplasmic diffusion coefficient, D mem (mm 2 /s) is the lateral diffusion coefficient along the membrane, k on,p (s À1 ) is the rate of photoconversion from dark to light state, k off,p (s À1 ) is the rate of thermal reversion from lit to dark state, k off, dark (s À1 ) is the rate of unbinding from the membrane in the dark state, k off,lit (s À1 ) is the rate of unbinding from the membrane in the lit state, k on,dark (M À1 s À1 ) is the rate of binding to the membrane in the dark state, k on,lit (M À1 s À1 ) is the rate of binding to the membrane in the lit state, and S max (molecules/mm 2 ) is the maximal surface density of binding sites. Importantly, all constants are empirically measured in this work ( Figure 2B), as detailed later.
Given the non-linearities in the PDEs, model initial conditions were expressed in absolute protein concentration units rather than relative fluorescence units, to allow for direct correlation to single-cell data and inter-cell comparisons. The conversion between mCherry fluorescence to BcLOV4 concentration was calibrated by established methods with intracellular access via (B) Binding constants and rates directly measured in live HEK single cells by cytoplasmic depletion fluorescence microscopy. Diffusion rates were determined by FRAP microscopy ( Figures S3 and S4). When needed for feasibility, upper/lower bounds were chosen as conservative values, as noted. Flavin photochemistry was measured by absorbance spectroscopy using FPLC-purified bacterial-expressed recombinant protein without fused mCherry (STAR Methods; l = 455 nm). The parameter set, along with fluorescence-derived absolute protein concentration ( Figure S1) and measured hardware PSF (STAR Methods), is sufficient to fit the global minimum of the proposed single-cell models.  (Cherkas et al., 2018). The cytosolic concentration was typically $0.5-3 mM protein (mean = 4.02 mM, 95% confidence interval [CI] = 3.44-4.82 mM, N = 107 cells). The mean of the calibrated single-cell imaging distribution mirrored the population average estimated directly from cell lysate (mean = 3.72 ± 0.79 mM, N = 4 plates) ( Figure S1). The requisite bulk-surface coupling cannot be accurately handled by built-in MATLAB routines or the standard versions of many multiphysics packages. Thus, we derived the weak form and constructed FEM matrices via FELICITY (Walker, 2018). A custom solver (see STAR Methods) was built for the purpose of computational efficiency. Run times ranged from 5 to 10 min per cell (per single core of a 3.7 GHz processor) for a typical experiment (duration = 200 s, Dt = 0.1 s, and maximum element size = 5 mm 3 ). Under these 3D mesh conditions, the error was <1% versus the outputs simulated with ultrafine mesh conditions of millisecond resolution on the 10-nm scale for a typical cell ( Figure S2). This throughput facilitated the simulation of every experimental cell for the high-resolution studies herein.

Measurement of biophysical constants
We experimentally derived all constants ( Figure 2B) necessary to validate the model and gauge its accuracy. The photoinduced reaction-diffusion model constants capture the flavin photochemistry, diffusion, and membrane binding kinetics. All parameters were measured in situ in HEK cells, except for photochemical constants measured by in vitro spectroscopy using purified recombinant protein that lacked the mCherry tag.
Diffusion constants were measured by fluorescence recovery after photobleaching (FRAP) microscopy ( Figure S3). The measured cytosolic diffusion (D cyto = 8.932 mm 2 /s) is consistent with our previous estimate (Glantz et al., 2018). The membrane lateral diffusion (D mem % 0.028 mm 2 /s) is slower than those of proteins bound to a freely mobile lipid or lipid-modified binding partner (0.1-1 mm 2 /s), but unsurprising for a protein that binds the membrane via distributed electrostatic interactions (Kenworthy et al., 2004;Hammond et al., 2009;Ambroggio et al., 2010;Ribrault et al., 2011;Yasui et al., 2014;Bendezú et al., 2015;Natwick and Collins, 2021). The low D mem indicates that BcLOV4 primarily spreads along the membrane not by lateral diffusion but rather unbinding and rebinding (or ''hopping' ' [Yasui et al., 2014]) events that cause spreading via cytosolic diffusion.
Membrane binding affinity and kinetics of BcLOV4 to the inner leaflet were measured in cellulo. The lit-state affinity was K D,lit = 1.19 mM (N = 23, Figures S4A-S4D), which is in line with similarly in situ-measured membrane affinities for other anionic lipid binding domains (Smith et al., 2015). The affinity in the dark-adapted state was K D,dark > 20 mM ( Figures S4E-S4G), which is reported as a minimum bound due to the lack of observable dark-state membrane localization such that k on,dark could not be determined by our bulk (versus single-molecule) technique in live cells. The maximal surface density (S max ) of BcLOV4, which governs its achievable dynamic range and the (non-)linearity of its membrane recruitment, was measured by monitoring maximal achievable cytosolic depletion; the membrane density is calculated from the protein that leaves the cytoplasm given the known surface area-to-volume ratio from the 3D-reconstructed confocal stacks. The measured value of $6,000 molecules/ mm 2 (N = 83) is similar to ones measured in situ for overexpressed proteins that bind endogenous membrane lipids (Smith et al., 2015).
PSF convolution for comparing models with real-world data Finally, we convolved the model output with the optical hardware PSF to improve direct comparison of model outputs to microscopy findings. Likewise, the optical input of the framework also factors in the diffraction-limited excitation beam profile to account for off-focal plane excitation, which can impact the signaling induced by cytosolically diffusible proteins.
Owing to diffractive effects, observed membrane fluorescence is inherently a partial function of cytosolic protein concentration, where the $7-nm-thick membrane accounts for just a fraction of the $100-400 nm axial depth-of-field of optical microscopes at high magnification. Thus, measurements of dynamic localization that equate membrane density with membrane fluorescence will vary widely with the hardware system. (Figures S5A and S5B). However, fluorescence quantification of the larger cytoplasm can be less susceptible to these confounds because its fluorescence is not as affected by diffractive effects or subtle morphological changes ( Figure S5C). Hence from here on, measurements of protein depletion and recovery are done by tracking cytoplasmic fluorescence.

Single-cell model performance
We compared model outputs with results from spinning-disk confocal stimulation and imaging of BcLOV4 ( Figure 3). The model was able to mirror the salient spatiotemporal behavior of BcLOV4, including undulations in cytosolic protein level from repetitive unbinding/rebinding cycles in response to sparsely pulsatile stimulation ( Figures 3A and 3B). Deconvolution of confocal z-stack data before geometry generation did not significantly impact resultant cell meshes or cell-wide mean-squared error, likely due to the excellent axial sectioning capabilities of confocal microscopy.
Like any other computational approach, the model's accuracy depends on the granularity of the initial mesh conditions. For example, the standard output does not account for shared membranes with neighboring cells or for lysosomal voids that appear as fluorescent puncta upon post-induction protein binding ( Figures 3A and 3B). Accounting for lysosomal binding had limited impact on BcLOV4 membrane recruitment accuracy at normal biogenesis levels, but slowed the simulation by $10-fold, scaling by O(n 2 ) the number of lysosomes ( Figure S6).
The 3D non-linear modeling approach was more accurate than 3D linear and 2D (non-)linear modeling (Figures 3C and S2D-S2F; see STAR Methods for calculation of error). Note that this comparison represents a best-case scenario for the 3D linear and 2D (non-)linear models, since we initialized them with the absolute protein concentrations rather than with relative or assumed concentrations (Valon et al., 2015;Natwick and Collins, 2021). These data suggest that both the geometric aspects and the non-linearity contribute to model performance. They also suggest that linear models of membrane recruitment that treat the membrane as an infinite binding sink ( Figure 3D), while valid at low cytosolic concentrations, are less appropriate for Cell Reports Methods 2, 100245, July 18, 2022 5 Article ll OPEN ACCESS overexpressed proteins in engineered cells. This influence can specifically be seen by the poor performance of linear models for high-concentration cells ( Figure S2F).

Intrinsic and extrinsic determinants of spatial confinement
We next used our single-cell approach to explore how optical stimulation paradigm and hardware affect spatiotemporal resolution, specifically the spatial confinement of membrane recruitment, a motivating metric in optogenetics. In brief, a small 1.5 3 0.5 mm stimulatory region of interest (ROI) was patterned by laser-scanning confocal microscopy (LSCM), and the resultant fluorescence membrane profile over 1 min was fit with a Gaussian to estimate the membrane spread ( Figures 4A-4C). The spatial spread primarily occurs during the initial few seconds post excitation when the protein rapidly diffuses in the cytoplasm (Figures 4C and 4D,N = 14 cells). Upon binding, the spreading occurs by lateral diffusion and by rebinding events subject to cytosolic diffusion. While great attention in optogenetics has been given to minimizing lateral diffusion to reduce spreading (Valon et al., 2015;Van Geel et al., 2018), its contributions to spatial resolution are likely outweighed by cytosolic diffusion when the photoactivated protein is cytosolic, as is commonly the case.
To systematically explore the determinants of spatial confinement, we simulated both the molecular-level biophysical interactions and the optical hardware used for induction and visualization ( Figures 5A and 5B). We predicted the sensitivity of spatial confinement to these two parameter classes for a typical protein concentration of 1 mM ( Figures 5C and 5D). For the biophysical parameter sensitivity analysis, we used the cellular geometry from the confocal experiment of Figure 4 and a range of biophysical constants informed by experiments here or elsewhere (Bartelt et al., 2018;Natwick and Collins, 2021). For hardware analysis, we used a geometry generated by hemi-ellipsoid projection of an idealized cell (see STAR Methods). cell, which lacks perceivable membrane unbinding or cytosolic repletion during each 1 s period. Model recapitulates nuclear void and subcellular distribution in the cytosol and plasma membrane but does not account for neighboring cells or for lysosome binding that greatly increases computational time with minimal accuracy benefit ( Figure S6). Scale bar, 5 mm. (C) Mean-squared error (MSE) of the described 3D non-linear model (of panels A and B) versus a 3D linear model, 2D non-linear model, and 2D linear model. Pooled dataset contains cells stimulated at 4 = 0.67%-10% (N = 17). Linear models performed worse for high protein concentrations ( Figures S2D-S2F). Paired Wilcoxon signed rank test: *p < 0.05, **p < 0.01. (D) Schematized non-linear modeling of dynamic membrane recruitment. The system is linear at low membrane binding site occupancy or low cytosolic BcLOV4. High concentrations typical of overexpressed inducible/optogenetic signaling systems require non-linear models due to high fractional occupancy of membrane binding sites upon photoactivation.
6 Cell Reports Methods 2, 100245, July 18, 2022 Article ll OPEN ACCESS Spatial resolution was most profoundly improved by parameters (lower D cyto , higher k on,lit , lower k off,lit , and higher k off,p ) that directly decrease the average distance traveled by a cytosolically diffusing lit-state protein by slowing it, improving binding, decreasing rebinding, or decreasing the lit-state lifetime (Figure 5C). An ideal D mem = 0 can only improve the confinement to the spatial resolution limit set by the other parameters . However, D mem is indeed limiting for the less common heterodimerization configuration where the photosensory protein is the membrane-localized component (Toettcher et al., 2011;Natwick and Collins, 2021). Of note, increasing S max improves resolution by reducing diffusion times and distances; it also enhances the maximum possible inducible signaling (Natwick and Collins, 2021).
The analysis of how optical hardware excitation volume impacts recruitment spatiotemporal dynamics found that a laserscanning confocal system (1P) induces greater signaling than more axially precise total internal reflection fluorescence (TIRF) and 2P, by exciting a proportionally larger volume for similar membrane surface areas ( Figure 5D). When only one optical pulse is delivered, TIRF and 2P potentially become less focal, with increasing SD as membrane-bound protein continues to spread after the excitation volume depletes of photoactivated protein (greater than tens of seconds); conversely in 1P, cytosolic photoexcited proteins continue to populate the illuminated surface area. If confinement is defined as the profile width at half the experimental maximum (instead of instantaneous relative maximum), TIRF and 2P are more focal because the protein den-sities are rapidly clipped by the quantification threshold. Thus, the primary spatiotemporal benefit of axial confined TIRF and 2P is an improved impulse response of induced signaling.
In silico resolution of inter-experiment variability Given its accuracy in predicting cell-wide dynamics, the framework provides a powerful means to improving robustness in single-cell studies by bridging inter-experiment variability in silico. To this end, we explored differences between subcellular digital micromirror device (DMD)-stimulation experiments on widefield versus confocal microscopes. We observe a phenomenon on widefield microscopes whereby cytoplasmic fluorescence appears enhanced within a patterned stimulation ROI, as if it continuously recruits protein from distal unstimulated regions ( Figure 6A). The paradoxical enhancement is the opposite of the expected rapid cytoplasmic depletion from the ROI. It is not observed when we use a confocal microscope to image a similarly DMD-stimulated cell ( Figure 6B). This discrepancy suggests that the enhancement is due to the poor axial sectioning capability of a widefield microscope ( Figure 6C). Post hoc analysis found that the hardware-dependent imaging artifact stems from the combination of (1) PSF-limited photonic integration of membrane-bound protein within the ROI and (2) simultaneous photobleaching throughout the cell ( Figures 6D-6F).
Next, we explored how cell geometry affects the observed kinetic constant for membrane association (t on ), a common metric for comparing optogenetic recruitment tools (Buckley et al., 2016;Zimmerman et al., 2016). While it is known that cell size (C) SD of the Gaussian profile (solid black line ± 95% CI, N = 14) from the heatmaps, a proxy for the degree of spatial confinement. (i) Experiment and (ii) model prediction. The SD initially rapidly increases because of cytosolic diffusion-limited association. BcLOV4 is then largely immobile due to slow lateral diffusion and membrane affinity. (D) Biophysical processes (i.e., excluding hardware contributions) that govern spatial resolution of optically inducible recruitment, pre-/post-binding to the membrane. See Figure 5 for sensitivity analysis.
Cell Reports Methods 2, 100245, July 18, 2022 7 Article ll OPEN ACCESS and geometry can impact signaling (McBeath et al., 2004;Meyers et al., 2006;Neves et al., 2008;Hara and Kimura, 2009), these factors are often overlooked in optogenetics. Spe-cifically, the plasma membrane surface area to cytoplasmic volume ratio (SA/V) critically impacts the binding site availability S max relative to total number of diffusible molecules. Increasing (2) unbinding frequency and distance traveled before rebinding, and (3) distance traveled by lateral diffusion. (C) Sensitivity analysis of spatial confinement to the intrinsic biophysical parameters. Values span those derived here for BcLOV4 and elsewhere for heterodimerization systems. SD of the protein distribution along a modeled membrane matching the experimental conditions of Figure 4. (i) Binding site availability (S max ), (ii) k on,lit , (iii) k off,lit , (iv) k off,p of the chromophore photocycle, (v) k on,dark , (vi) k off, dark , (vii) cytosolic diffusivity, or (viii) lateral diffusivity along the membrane. The role of lateral diffusion is limited for large excitation volumes and/or when binding kinetics permit extensive cytosolic diffusion when rebinding. (D) Effect of stimulation method. Systems with larger excitation volumes drive more potent signaling but reduce spatial precision. (i) Simulated excitation of 2 3 2 mm region at the bottom of a model cell (excitation duration = 100 ms, cytoplasmic [BcLOV4] = 1 mM) by different stimulation methods (10 W/cm 2 irradiance at focal plane). Highlighted volumes show the 1 W/cm 2 isosurfaces of the excitation volume. (ii) Peak membrane-bound protein recruitment, baselined to subtract pre-bound protein in the dark state. The rapid drop-off in peak recruitment for axially confined stimulation (2P, TIRF) due to low levels of activated cytosolic BcLOV4 results in decreased induced signaling. (iii) Spatial resolution quantified by SD for a Gaussian distribution fit (1P) or the full width at half maximum /2sqrt(2ln(2)) (2P, TIRF) of the membrane profile at each time point. Axially confined methods do not outperform classic 1P excitation for this measure overall, but better retain focality during initial recruitment (t < 5-10 s). (iv) Spatial resolution quantified by width at half of initial/absolute maximum. Inset schematizes the steep drop-off at longer time points by this metric from threshold clipping. We also demonstrated that a cell-unique model experimentally validated in one experiment can be used to predict the outcome of another of entirely different hardware, stimulation paradigm, and analysis ( Figure 7C). We compared t on determined by (1) fitting cytosolic depletion from confocal micrographs of pulsatile-stimulated cells (i.e., method used here) versus (2) colocalization correlation analysis with a membrane marker from widefield images of continuously stimulated cells (our previous method [Glantz et al., 2018]). Notably, the marker here is purely virtual, not coexpressed as before. In silico transposition recapitulated a salient interpretive confound that t on derived by cytoplasmic depletion ($5 s) is systematically slower than by marker colocalization ($1.5 s). The latter perceived acceleration is a consequence of the mathematical correlation ( Figure 7D). Because the marker has background cytoplasmic fluorescence, the acceleration reflects the synergistic increase in correlation along both the cytosolic and membrane segments and a faster approach to mathematical saturation. Such predictive simulations show how computational approaches can improve interpretive robustness by permitting results from one lab or experiment to be compared with those emulating another.

Implementation in other computational environments
While the mesh generator and solver were custom-built for computational efficiency, the single-cell framework can be implemented in other software. It is modular in that the (1) 3D cellular mesh, (2) time-varying 3D-patterned distribution of optical excitation, (3) reaction-diffusion model (or other relevant model), and (4) PSF convolution can all be treated as separate components imported and exported between computational environments. Thus, we implemented the approach in VCell (Figure S7), a general environment for simulating cell signaling (Moraru et al., 2008) that can natively mesh a 3D cell (based on the finite volume method, with a uniform orthogonal mesh) and account for PSF convolution. We integrated the calculation of spatiotemporally complex 3D distributions of optical excitation (e.g., laser raster) into VCell by introducing the optical input distribution computed within our framework as a specific ''light'' element, taking inspiration from work elsewhere in modeling FRAP experiments (K€ uhn et al., 2011) ( Figure S7A).
We successfully replicated our analyses in an end-to-end fashion with closely approximated results ( Figures S7B-S7D, related to Figure 5D). The computational runtime increased in VCell versus the fully custom solution ($10-fold slowdown to simulate a TIRF experiment). Since its inputs are continuous functions, VCell recomputes the intensity of excitation light at each meshpoint at every timestep, even when it is unchanged. Conversely, the custom FEM can utilize coarser mesh sizes to reconstruct the topologies and store the optical input elements in local memory. Nevertheless, the successful replication demonstrates that one can employ the modeling tools of choice and those best suited for the single-cell modeling task at hand.

DISCUSSION
Predictions by our single-cell framework agreed with experimental single-cell data. Many features contribute to its utility in studying peripheral membrane protein dynamics: 3D nonlinear dynamics without concentration assumptions, volumetric reconstructions of single cells, optical beam profiling and PSF convolution to simulate hardware-specific effects, and validation by a global optimal solution for the set of PDEs (i.e., all the fundamental constants and absolute concentrations were measured).
While we successfully simulated anticipated outcomes and visual outputs of single-cell experiments, the quality of the simulated micrograph depended (1) on the mesh fidelity of the initial volumetric reconstruction as already described, and (2) on the underlying biophysics accounted for in the model. For example, the model does not recapitulate BcLOV4 recruitment ''hot spots'' (e.g., punctate membrane fluorescence in Figures 1  and 3A), which likely form by de novo aggregation on the membrane given that BcLOV4 is monodisperse in the cytosol (based on identical FRAP-measured and theoretical cytosolic diffusivities). While we cannot exclude a role for curvature in this phenomenon, our current understanding from in vitro studies (of lipid headgroup selectivity) of the BcLOV4 membrane interaction is that it is more governed by membrane anionic charge density than by curvature (Glantz et al., 2018).

Optogenetic membrane recruitment-based signaling implications
Single-component BcLOV4 recruitment reaches a fractional occupancy of [S]/S max $ 30% for typical cellular concentrations, similar to the occupancy reached by the iLID system (Natwick and Collins, 2021). Since the total induced signaling magnitude depends on the cytosolic concentration of the payload, a Article ll OPEN ACCESS single-component system whose surface-bound partner is in natural excess can decrease the metabolic load by $80% over an optimized heterodimerization system, i.e., one expressed with a 4-fold excess of the surface component (Natwick and Collins, 2021). Likewise, since the BcLOV4 S max ($6,000 molecules/mm 2 ) exceeds that of lipidated peripheral membrane proteins (500-1,500 molecules/mm 2 [Del Piccolo and Hristova, 2017;Angert et al., 2020]), BcLOV4 likely has a greater achievable membrane density than a heterodimerizer; consequently, it may drive greater signaling induction over a broader concentration range of linear binding kinetics.
It was unanticipated that when the photosensory domain is cytosolic, its cytosolic diffusion has a much larger role than lateral diffusion in the spatial confinement, which becomes especially important under binding-limiting cases. Also unexpectedly, we found that axially confined stimulation improved the impulse response of induced signaling more than the spatial resolution. However, in systems with a membrane-localized photosensor, the confinement should be determined primarily by lateral diffu-sion and is largely independent of excitation volume. Thus, a system with a membrane-localized photosensor is likely more ideal for maximizing spatial confinement; however, this configuration comes at the expense of increased metabolic load and/or decreased signaling magnitude for the fractional occupancy reasons described above.
Our studies elucidated the biophysical determinants of performance in optogenetic membrane recruitment and emphasized how in silico approaches can enhance interpretive robustness in single-cell signaling studies. The validated framework provides a foundation for the principled design of improved tools and complex stimulation paradigms for tight and potentially closed-loop spatiotemporal control of cell signaling. Considerable progress in optogenetics has been made in model-guided signaling gradient formation (de Beco et al., 2018;Izquierdo et al., 2018), and improved strategies could further advance the field in the study of oscillatory waves and more complex pattern formation (Chiou et al., 2018  BcLOV4 refers to BcLOV4-mCherry. Kinetics measured in one experiment can be reasonably extrapolated to data generated in another (e.g., different cell or experimental condition). (A) t on varies across a physiological range of cell surface area-to-volume ratios (SA/V) that determine membrane binding site availability relative to total intracellular BcLOV4 (0.1 s pulse, l = 405 nm, duty cycle 4 = 10%, 12.24 W/cm 2 ). SA/V range of $0.15-1.2 mm À1 spans geometries from a large 20-mm-radius cell with negligible nuclear fractional volume to a small 5-mm-radius cell with sizable 4-mm-radius nucleus. (B) Log-order acceleration of t on with increasing SA/V as the main contributor to intercellular differences in recruitment kinetics. Experimental cells and corresponding FEM (experiment: black, 95% CI; simulation: red) and calculated t on (blue band, in idealized spherical cells) across SA/V and stimulation duty ratios (as in A, 4 = 0.67%-10%). (C) In silico ''transposition'' of a cell between experiments. (i) Validated cell-unique FEM from a pulsatile stimulation experiment with t on quantified by cytosolic depletion is simulated as if it was in a different experiment of continuous stimulation and quantified by colocalization with a virtually introduced membrane marker. (ii) Simulated transposition of one cell with virtual GFP-CAAX (experiment-derived marker background fluorescence of $10% relative to its membrane fluorescence; error = 95% CI; $0.5 mM protein to match previous work [Glantz et al., 2018]). Colocalization correlation analysis along a line profile transecting the membrane results in faster perceived recruitment dynamics than by cytoplasmic depletion. Predictions on markerless experimental cells using virtual markers here agree with previous work (inset table). (D) Post hoc analysis reconciles data quantified by the two methods. (i) Schematized GFP-CAAX marker with constant membrane/cytosol ratio. (ii) ToolTmarker correlation initially improves in both the membrane and cytosol to synergistically accelerate perceived colocalization increase. The correlation improvement subsequently slows when membrane colocalization increase is counterbalanced by cytosolic colocalization decrease (when cytosolic BcLOV4 depletes beyond the cytosolic marker level).
Cell Reports Methods 2, 100245, July 18, 2022 11 Article ll FEM in cell signaling FEM is widely used to study mechanics and transport, including in cell motility (Neilson et al., 2010;Elliott et al., 2012) and deformation biomechanics using non-stationary surface elements (Boccaccio et al., 2011;Hughes et al., 2018). Its application for spatially modeling intracellular signaling is less common but more theoretically advanced for describing a coupled bulksurface system of membrane recruitment than alternative approaches (Madzvamuse and Chung, 2016;Cusseddu et al., 2019;Paquin-Lefebvre et al., 2020). Beyond equation-based applications, rules-driven or agent-based modeling frameworks such as PhysiBoSS (Letort et al., 2019) are particularly useful when spanning multiple scales, such as incorporating intracellular signaling cascades with cell-cell interactions (Bonabeau, 2002;Gorochowski, 2016). Most existing multi-scale agentbased frameworks of cell signaling still approximate cells as ideal spheres. We envision that single-cell FEM and excitation light beam computation techniques described herein are extendable to rules-driven models and would benefit them by facilitating tolerance to experimental imaging conditions and cell heterogeneity.
In summary, we described a conceptual framework that can improve computational models of single-cell signaling and microscopy experiments. We applied it to study BcLOV4 spatiotemporal dynamics to identify key performance determinants in optogenetic membrane recruitment-based signaling. We demonstrated how the approach can be used for improved robustness when interpreting single-cell data. By accounting for both the intrinsic biophysical contributions of the optogenetic tool and the extrinsic contributions of the optical system used for stimulation and imaging, this approach can guide forward the design of improved optogenetic tools and experiments for dissecting peripheral membrane protein signaling dynamics.

Limitations of the study
The current framework was designed for single-machine computation. Future work includes expanding its throughput via distributed computation functionality. While VCell models can be distributed to run on its cloud server, the current study was not compatible with this online server owing to the presence of spatiotemporally varying inputs. Specifically, VCell does not natively provide a simple interface through which to model time-and space-varying inputs at this time. The described local VCell implementation introduced this functionality without access to the underlying codebase; this approach required single-machine computation. Accordingly, the described runtime performance assessments represent hardware-matched headto-head comparisons. Despite these limitations, the described framework is more broadly envisioned as an adjunct to existing tools and is intended to streamline the analysis of complex single-cell signaling studies.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

ACKNOWLEDGMENTS
We thank Christopher Fang-Yen for helpful discussion, Andrea Stoudt for technical assistance, and Lukasz Bugaj for confocal microscope access.

DECLARATION OF INTERESTS
The authors declare no competing interests.

EXPERIMENTAL MODEL AND SUBJECT DETAILS
Cell lines HEK293T (ATCC, CRL-3216) cells were cultured in D10 media composed of Dulbecco's Modified Eagle's Medium with GlutaMAX (Invitrogen 10566016), supplemented with Penicillin-Streptomycin at 100 U/mL and 10% heat-inactivated fetal bovine serum (FBS). Cells were maintained in a 5% CO 2 tissue-culture treated dishes in a water-jacketed incubator (Thermo/Forma 3110) at 37 C.

METHOD DETAILS
For the following sections, unless otherwise specified: (i) water was Milli-Q water (ddH2O, 18.2 MU $ cm), (ii) oligonucleotides were synthesized by IDT, (iii) genetic constructs were verified by Sanger sequencing (Genewiz), (iv) all experiments referring to BcLOV4 specifically refer to the BcLOV4-mCherry fusion construct, unless specifically stated otherwise. Native sequence and mammalian codon-optimized BcLOV4 plasmids are publicly available via Addgene.

Genetic constructs, protein expression and transfection
Recombinant BcLOV4 production BcLOV4 was bacterially produced and purified by FPLC (fast protein liquid chromatography) as described previously (Glantz et al., 2018). Protein was stored at 4 C and used within 2 days of purification. Photocycle kinetics and quantum efficiency BcLOV4 (without fused mCherry) absorbance scans were used to determine photocycle kinetics by monitoring absorbance at 450 nm (A450) as previously described (Glantz et al., 2018). To measure quantum efficiency of photoconversion (4Þ and k on,p , 15 s of baseline measurements were made at room temperature, and the A450 immediately was recorded after blue light stimulation where we can ignore lit-to-dark state thermal reversion because it occurs on timescales much longer than photoactivation; the fraction of activated BcLOV4 after 5 s of excitation will be 1 À e À 5kon;p : The constants k on,p and 4 are related to each other by the rate of photon absorption: where s ex is the absorption cross section in units of cm À2 , 4 p is the photon flux density in units of photons cm À2 s À1 , and l is the wavelength of excitation light. Given: where e is the molar extinction coefficient ($12,500 M À1 s À1 ) and N A is Avogadro's number, and: where P d is the power density of the excitation light, h is Planck's constant, and n is the frequency of excitation, we approximate this integral by assuming that our LED spectrum is narrow and centered at 455 nm: The fraction of BcLOV4 that is activated after 5 s is: which can be fit to the experimental data to extract 4 and then use the value to calculate k on;p . Dark-state reversion kinetics were measured by monitoring the absorbance at 450 nm (A450) post-stimulation for every 0.5 s for 2 min, and k off,p was fit as an exponential to reversion kinetics. Mammalian cell culture and transduction HEK293T (ATCC, CRL-3216) cells were cultured, plated, and transfected as previously described (Glantz et al., 2018). Cells were imaged 24 h after transfection.

Microscopy calibration and assays
Widefield microscopy and spatial patterned excitation Widefield imaging was performed on an automated epifluorescence microscope equipped with a custom-built spatial patterning system (Berlew et al., 2020) to deliver blue light (l = 455/20 nm, 12 mW/cm 2 , 100-1000 ms pulses at duty cycle = 0.8-10%). Excitation patterns were typically 25 mm-wide squares that illuminated $25-50% of selected cell area. Spinning disk confocal microscopy A Nikon Eclipse Ti2 microscope equipped with a CSU-W1 confocal scanner unit (Yokogawa) was used with a 60x/1.4 NA objective. BcLOV4 was stimulated at l = 405 nm at either 50% [12.24 W/cm 2 at the focal plane] or 100% [16.94 W/cm 2 at the focal plane] laser power. mCherry was excited at l = 552 nm and imaged using a Brightline mCherry-C-000 filter cube (562/40 excitation filter, 593 dichroic, and 641/75 emission filter). When necessary, spatially patterned excitation could be delivered via a Nikon DMD module (tunable irradiance between 10 and 70 mW/cm 2 at the focal plane). In general, excitation intensity was set at $12 mW/cm 2 to allow for accurate inter-comparison between widefield and confocal DMD excitation datasets. Excitation intensity measurements For widefield wholefield, widefield DMD excitation, and confocal DMD excitation, intensity at the focal plane was measured with a Thorlabs PM100D power meter. PyOptica (Grochowicz and Miler, 2020) was used to calculate the expected out-of-plane excitation or light intensity at different z-locations that would yield the observed DMD image formation at the focal plane. Because of the underfilling of the back-aperture of the objective, the excitation volume for DMD excitation is, at mm-distance scales, a vertical square column of light with identical irradiance profiles in each z-plane. This approximation was used for all DMD simulations.
For wholefield spinning disk confocal excitation, illumination intensity at the focal plane was measured by fitting photobleaching curves of <5 mm water-in-oil emulsion droplets filled with 100 mM riboflavin. To estimate out-of-plane excitation, we approximated the spinning disk unit as continuously scanning multiple diffraction-limited Gaussian beams across the e3 Cell Reports Methods 2, 100245, July 18, 2022 Article ll OPEN ACCESS field-of-view (neglecting any beam-shaping occurring in the unit). The intensity of a Gaussian beam can be modeled as (Dickson, 1970): where r = ffiffiffiffiffiffiffiffiffiffiffiffiffiffi x 2 + y 2 p and w is the beam radius and is defined as: where z R is the Rayleigh range: where n sample is the refractive index of the propagating medium and l is the wavelength of the light. Note that w 0 can be estimated from the point spread function (calculated as 109.8 nm). The time-averaged volumetric intensity distribution of delivered excitation light was then approximated as the convolution of I(r,z) with a 2D kernel representing the amount of time the spinning disk scanner spends at each point. For whole-field excitation, the average light intensity that out-of-focus planes of the cell are exposed to was found to be approximately equal to the light intensity at the focal plane and hence a uniform excitation light intensity throughout the cell volume was assumed in corresponding models.

Point spread function measurements
Spinning disk confocal PSFs were measured by taking z-stacks of 100 nm fluorescent beads (Degradex 2101C/2211). Widefield PSFs were measured using larger beads 100-390 nm fluorescent beads (Spherotech FP-0252-2/FCM-02556-2). PSFs were extracted from z-stacks using the Python Microscopy Environment (Baddeley et al., 2021) and subsequently de-convolved to account for non-negligible bead diameter. 3D Gaussian (spinning disk) or Richards & Wolf (widefield) PSF models were fit to the experimentally derived PSFs with the python PSF package (Gohlke and Holub, 2021) to limit noise for subsequent analysis.

In cellulo measurement of kinetic constants
The values for k on,lit , k off,lit , k on,dark , and k off, dark , were extracted from in situ HEK cell measurements. Given the fast cytosolic diffusion relative to the membrane binding rate of BcLOV4, binding was modeled as reaction-limited, and assuming isotropic diffusion, was simplified to a four-state zero-dimensional model: where L is a characteristic length-constant, defined as the volume-to-surface area ratio of the cell, which is necessary to maintain conservation of mass in the system. L is directly measurable via confocal imaging. For BcLOV4 membrane unbinding post-excitation when t [ 1 k off;p , u M ; u C z0 and the model can be simplified to: Due to the relatively small expected value of k on,dark , it is neglected to recover the expression: Based on simulations, the model was accurate to within 5% for t > 40 s post-completion of excitation. k off, dark was thus determined from a spinning-disk dataset of cytosolic fluorescence recovery after 5 s wholefield excitation of BcLOV4 at 16.94 W/cm 2 ( Figures S4E and S4F). To fit k on,lit and k off,lit , the 0D-model in Equation (12) was fit to a spinning disk dataset of dark-adapted BcLOV4 stimulated 100 ms every 1 s (10% duty cycle) at 16.94 W/cm 2 (Figures S4A-S4D). Due to relative model insensitivity to the value of k on,dark , which primarily affects dark-state affinity and has minimal effect on membrane binding or unbinding kinetics, it was better estimated as a bound. Given the fact that BcLOV4 at the membrane could not be visualized in the dark state regardless of cytosolic BcLOV4 concentration, it followed that membrane fluorescence (F mem ) Cell Reports Methods 2, 100245, July 18, 2022 e4 Article ll OPEN ACCESS was always less than or equal to cytosolic fluorescence (F cyto ), i.e., F mem % F cyto . F mem is approximately a linear function of the protein density at the membrane, which in turn is related to the dark-state BcLOV4 concentration by the Langmuir isotherm. Hence: A 1 , A 2 were determined using the giant unilamellar vesicles (GUV) calibration (see ''Surface density calibration''). Similarly, for F cyto : where B 1 , B 2 were measured with patch micropipette microinjection of Lucifer Yellow (see ''Concentration calibration''). Assuming F mem % F cyto , the lower bound on K D,dark of 20 mM was determined from Equations 15 and 16, and consequently the upper bound on k on,dark of 1125 M À1 s À1 , was used in the model as a ''worst-case'' value ( Figure S4G).

Surface density calibration
The conversion between BcLOV4 surface density and membrane fluorescence was measured based on techniques by others (Dezi et al., 2013). Briefly, phosphatidylcholine giant unilamellar vesicles (GUVs) with varying concentrations of Texas Red-1,2-dihexadecanoyl-sn-glycero-3-phosphoethanolamine, triethylammonium salt (Tx-DHPE, Biotium) between 0.05 and 1% mole/mole were generated using the water-in-oil emulsion transfer method (Nishimura et al., 2012). GUVs were imaged at their equatorial plane by spinning-disk confocal microscopy, and a calibration curve was constructed which related Tx-DHPE density to measured membrane fluorescence. GUVs were automatically identified in the image by the circle Hough transform and membrane fluorescence was extracted in MATLAB as the average signal measured along the circumference. The relative brightness of Texas Red versus mCherry under our imaging conditions converted this to a curve of BcLOV4 density versus membrane fluorescence. On the spinning disk confocal microscope used here, the measured membrane fluorescence is 1.2-fold the BcLOV4 surface density in molecules/mm 2 .
Concentration calibration Due to model nonlinearities, all inputs and outputs need to be expressed in absolute concentration units. The fluorescence-derived concentrations were measured based on previously used methods (Cherkas et al., 2018), wherein a known concentration of dye is delivered to cells by whole-cell patch-clamp ( Figure S1A). HEK cells' transmembrane potentials were kept stable using an Axopatch 200B amplifier and Digidata 1440 digitizer (Molecular Devices) at room temperature. Patch micropipettes were pulled to obtain a resistance of 3-10 MU (1.5 3 0.86 mm borosilicate glass, Sutter Instrument P-1000 pipette puller). Extracellular solution (Tyrode's solution) consisted of 125 mM NaCl, 2 mM KCl, 3 mM CaCl 2 , 1 mM MgCl 2 , 10 mM HEPES, 30 mM glucose, pH 7.3 (NaOH adjusted), 300 mOsm (sucrose adjusted). Intracellular solution consisted of 125 mM K-gluconate, 8 mM NaCl, 0.1 mM CaCl 2 , 0.6 mM MgCl 2 , 1 mM EGTA, 10 mM HEPES, 4 mM Mg-ATP, 0.4 mM Na-GTP, pH 7.3 (KOH adjusted), with 295-300 mOsm (sucrose adjusted). Intracellular solution was supplemented with Lucifer Yellow (LY) dye (5-25 mM) for perfusion. Cells with a leak current > 200 pA or access resistance >25 MU were discarded. Time-lapse fluorescent imaging was conducted at 6 Hz to monitor LY filling of the cell cytoplasm, which took $5 min to saturate. (N = 3 cells per 5-25 mM LY concentration). Known concentrations of BcLOV4-mCherry and LY were imaged in bulk solution, to derive the relative fluorescence ratio between these two fluorophores for the microscope. With this information, a linear function was constructed to approximate cellular BcLOV4 concentration (in mM) from fluorescence visualization intensity, which was given by F = 505½BcLOV4 + 125 for the widefield microscope and F = 452½BcLOV4 + 107 for the spinning-disk confocal microscope, given the imaging parameters used herein.
Based on these empirical relationships, the distribution of cellular BcLOV4 concentrations after transient transfection were determined from confocal data ( Figure S1B) to be a right-ward skewed distribution with a mean concentration of 4.02 mM. As secondary validation of this measurement, mCherry fluorescence in HEK cell lysate (SoluLyse-M, Amsbio) was determined from calibrated fluorescence of FPLC-purified recombinant BcLOV4-mCherry, after correcting for cell confluency, transfection efficiency, and cell volume/lysate volume; the confirmatory calculated mean concentration was 3.72 ± 0.79 mM (N = 4). FRAP measurements FRAP was performed using a Leica TCS SP8 laser scanning confocal microscope ( Figure S3). For lateral membrane diffusion, BcLOV4 was initially recruited to the membrane with a 100 ms pulse of 405 nm light. A small rectangular ROI (average size $1.5 3 0.5 mm, with the long axis parallel to the membrane) on the membrane was photobleached using a high intensity 561 nm laser at 100% power. The ROI was then imaged at $60 Hz until its average fluorescence stabilized. When fit to a 1D reaction-diffusion model that accounted for membrane unbinding/rebinding, the membrane unbinding/rebinding could account for most of the observed fluorescence recovery, resulting in D mem values below the detection limit. Therefore, the effective upper bound on D mem was determined by fitting the data to a canonical model without the reaction components (Ellenberg et al., 1997): where w is the width of the bleaching rectangle, and I is the intensity within the ROI. For measurement of cytosolic diffusion, a 2 mm diameter circular ROI near the center of the cell was photobleached using a high intensity 561 nm laser at 100% power ($250 ms bleach time). Fluorescence recovery was measured in the Leica ''fly mode'' to capture the initial few ms of the recovery curve and tracked at $600 Hz until the signal stabilized. The half-time of recovery was extracted from the data through fitting of an exponential recovery curve, to calculate D cyto using a previously reported conversion for circular ROIs (Kang et al., 2012). Spatial confinement Spatial confinement experiments were conducted using a Leica TCS SP8 laser scanning confocal microscope. A dark-adapted cell was imaged at $5 Hz for 4 s, after which a small perimembrane ROI ($1.5 x $0.5 mm) was photoexcited with a $200 ms pulse of 405 nm laser light. The cell was subsequently imaged at $5 Hz until the fluorescence in the ROI returned to baseline. The ROI and membrane profile fluorescence kinetics were extracted using a custom-written ImageJ macro. Results were compared to model predictions which utilized the experimental cell geometry. Gaussian curve fitting was conducted in MATLAB.

Simulations
Mesh initialization and quality control For all simulations, unless otherwise specified, all cytosolic mesh nodes were initialized with the cellular cytosolic concentration calculated from the input fluorescence micrograph. Concentration was assumed to be uniform in the cytoplasm because z-stacks analysis did not show variability that was not primarily attributable to organelle voids, cell geometry, or diffraction. Membrane mesh nodes were initialized as S max K D;dark =ðK D;dark + ½BcLOV4 cyto Þ. The volume of mesh elements was <5 mm 3 . Mesh size was selected so that, on average, a <1% error relative to super-fine mesh ($10 À4 mm 3 ) was maintained ( Figures S2A and S2B). Timesteps varied from 0.01-0.1 s and were selected so that a <1% error relative to a timestep of 0.001 s was achieved ( Figure S2C). As secondary validation, in all cases, it was verified that the resultant solution did not change with a finer mesh or smaller timestep. Spatial discretization using FEM The approach is based on bulk-surface reaction-diffusion systems (Madzvamuse and Chung, 2016;Cusseddu et al., 2019;Paquin-Lefebvre et al., 2020). As shorthand, S = S max À u M À v M : The weak form of this system defined by Equations (1) and (2) is (letting 4; j be our test functions on U and dU, respectively):  In matrix form, we can express this as: )

Time discretization and integration
For time integration, the fractional theta method was used, which is an implicit method with the A-stability of implicit Euler and second order temporal convergence of Crank-Nicholson (Chrispell et al., 2007;Madzvamuse and Chung, 2014). In this method, the desired time interval is first split into equal length steps, henceforth termed t. For each time step, 3 sub-steps are performed. The equations to be solved for each time step are outlined below, where q = 1 À 1 = ffiffiffi 2 p .
The first sub-step solves the following linear equations:

Cell segmentation
Cell nuclei and plasma membrane contours were detected using the MATLAB active contours algorithm (Chan and Vese, 2001). Briefly, images were initially manually cropped to include only individual cells. Nuclear locations were seeded by drawing a small rectangle in the center of the nuclear void and the nuclear contour was evolved using the rectangle as a starting mask. The plasma membrane contour was evolved by using the entire image as a starting mask. For confocal z stack segmentation, the nuclear/cytoplasmic contours from the previous image in the stack were used as the seeds for the following image. For images where cells were touching ($80% of widefield images and $50% of confocal images), it was necessary to manually crop-out the cells that were not of interest and seed the contours with a hand-drawn contour in ImageJ.

Generation of volumetric cell mesh
When taking confocal data 100 nm-step z-stacks of the cell volume were acquired before initiating any experiments to reconstruct a volumetric representation of the cell by interpolating between the segmented 2D frames of the z stack. After the cell volume was constructed, it was converted to a surface triangulation. The surface triangulation was then meshed via MATLAB's built-in mesh generation algorithm.
In the case of widefield images, which lacked reliable z-stacks because of poor axial resolution, the volumetric representation was constructed by a hemi-ellipsoid projection of the nucleus and cytoplasm from a single 2D frame ( Figure 1D). Briefly, the heights for the nucleus and cytoplasm were set as 2 and 3 mm, respectively. The (x,y) position of the nucleus or cytoplasm ellipsoid zenith and nadir was selected by taking the center-of-mass of a binary mask of the nucleus/cytoplasm. To account for the flat portion of the cell touching the cover glass, the lower portion of the projection was constructed as a hemi-barrel shape. This hemi-barrel was generated by sectioning, with a plane perpendicular to the z axis, a hemi-ellipsoid of height equal to twice the desired height of the hemi-barrel. After the cell volume was constructed, it was converted to a surface triangulation, which was then meshed via MATLAB's built-in mesh generation algorithm. Sensitivity to biophysical parameters Upper and lower bounds for D cyto , D mem , k off,p , k on,lit , k off,lit , k on,dark , k off, dark , and S max were chosen based on biophysical feasibility and by comparison to values derived for heterodimerization systems. Parameter values were then selected that spanned the space between bounds, yielding 31 total parameter sets where only one parameter at a time differed from BcLOV4's measured biophysical parameters. Using the same cellular geometry and excitation ROI from Figure 4 and an initial cytosolic protein concentration of 1 mM, the model was then run on each of the 31 parameter sets and analyzed as described in ''Spatial confinement.'' The parameter k on,p was excluded from the analysis as it does not directly affect spatial confinement. Sensitivity to optical hardware diffraction To compare the spatial resolution achievable by laser-scanning confocal, two-photon, and TIRF stimulation, a common cell volume was generated by hemi-ellipsoid projection (nuclear height 4 mm, cytoplasmic height 6 mm) of a 2D cell area, where the nucleus was a circle with r = 15 mm and the cytoplasm was a circle with r = 25 mm. Note that the very large cell ensured no artifacts were introduced by protein diffusing all the way around the cell while still accounting for computational memory constraints. Since a very fine mesh was required to ensure sufficient resolution to finely differentiate the excitation volume, the cell volume was meshed to a maximum element size of <0.025 mm 3 (typically 0.002-0.005 mm 3 ), which yielded >10 6 elements per cell and utilized $16 Gb of RAM.
The simulated excitation (l = 450 nm, 10 W/cm 2 at the focal plane) region was a 2 mm 3 2 mm square centered on the bottom flat portion of the cell volume to allow for reasonably similar conditions for the different hardware. To account for excitation beam rastering within the illumination ROI, one illumination spot for each 100 nm 3 100 nm square of the illumination ROI was assumed, for 400 total raster points. The rastering was simulated as being instantaneous since the scan speed is fast relative to the time-course of the experiment.
For each of the microscopy methods, the simulated axial extent and the diffraction-limited nature of the optical beam path were also modeled. TIRF was modeled as using a 63x/1.4 NA objective with an angle of incidence of q incidence = 79 and with refractive indices n sample = 1.33, n coverglass = 1.52; then the z-dependence of light intensity was modeled as (Fish, 2009): where d is given by: where l is the wavelength of excitation light. Equations 9, 10, and 11 were used to model the xy-dependence of the beam, noting that w 0 , or the waist radius, can be approximated as (Young et al., 2015): For the case of one-photon stimulation, specifically laser scanning confocal, Equations 9, 10, and 11 were applied to determine the axial dependence, except with the beam waist determined from the PSF as 109.8 nm. For the case of two-photon e11 Cell Reports Methods 2, 100245, July 18, 2022 Article ll OPEN ACCESS stimulation, the axial dependence modeling proceeded similarly, except the intensity was given as the square of Equation 9 and l 2P = 2l 1P = 900 nm.
For each of the microscopy methods, we modeled a 100 ms stimulation and simulated 60 s post-stimulation. Unlike for laser-scanning confocal stimulation, the membrane profiles generated via TIRF and two-photon stimulation were not Gaussian; hence, their respective spatial spreads were calculated as the full width at half maximum (FWHM) and then converted to an effective standard deviation via FWHM = 2 ffiffiffiffiffiffiffiffiffi 2ln2 p STD, which holds for a Gaussian distribution, to allow comparison to Figure 4 and the biophysical parameter sensitivity analysis. The FWHM relative to both the instantaneous maximum and the global maximum were calculated to highlight the impulse response modulatory nature of TIRF and two-photon microscopy. Lysosome binding Datasets were collected via spinning disk confocal stimulation and imaging of cells co-expressing transiently transfected BcLOV4 and the LAMP1-miRFP670nano (Addgene plasmid #127435) lysosomal marker. 5 mM biliverdin was added to the imaging medium immediately before visualization to increase miRFP670nano fluorescence; non-specific partitioning of the cofactor into the plasma membrane was subtracted during analysis.
Models that accounted for the presence of lysosomes required the segmentation of lysosomal position from collected z-stacks (150 nm steps) of LAMP1-miRFP670nano fluorescence using a custom MATLAB script. In brief, the initial z-stacks were first deconvolved with the microscope PSF to reduce the impact of diffraction on downstream analysis. The xy-center of each lysosome was identified by sequential analysis of each 2D micrograph in the z-stack and by fitting each bright lysosomal puncta with a 2D Gaussian distribution (where the center of each distribution was the center of the corresponding lysosome). Each lysosome was visible in several z-stack frames despite the de-convolution step, and thus, its z-center was taken as the middle frame from the series of frames it appeared in. For practicality, the radius of all lysosomes was taken to be 340 nm (Woldemichael and Rosania, 2017). For theoretical simulations of the impact of lysosome volume density on the binding kinetics of BcLOV4, a mesh from an actual cell was randomly initialized with a given volume density of lysosomes and meshed as described previously. Paradoxical fluorescence enhancement In cases where BcLOV4 was imaged by a widefield microscope, an accurate 3D cell volume was absent and required volumetric extrapolation. To explore model accuracy for these conditions, videos of DMD-excited cells were collected for 1-10% duty cycles of stimulation (100 ms every 10 s, 250 ms every 10 s, 1000 ms every 10 s, or 1000 ms every 15 s). For each stimulated cell, a volumetric mesh was constructed by 3D extrapolated projection (as described in ''Generation of volumetric cell mesh''), and 100 s of BcLOV4 spatiotemporal dynamics were simulated. To correct for photobleaching, biexponentials were fit for t > 150 s, i.e., when cellular concentration was at steady-state. To better visualize the recapitulation of paradoxical fluorescence enhancement within the stimulation ROI, the convolved model prediction was decomposed into its cytosolic and membrane components, since convolution is a linear operation.
To verify that the observed cytosolic fluorescence enhancement was truly an artifact of the diffractive properties of a widefield microscope, a similar experiment was conducted with DMD excitation and spinning disk confocal imaging. Specifically, data for DMDexcited cells was collected for 1-10% duty cycles of stimulation. In all cases, fluorescence enhancement within the excitation region was absent. Cell morphology influence on association kinetics Initially, meshes for spherical cells (r = 5-20 mm) with variable sized spherical nuclei (r = 1-10 mm) were generated. The resultant cell surface area-to-cytoplasmic volume ratios (SA/V) ranged from 0.15-1.2 mm -1 . Simulations were then run on these geometries for 1 mM initial cytosolic BcLOV4 for duty cycles ranging from 0.67-10%, an excitation duration of 0.1 s with 405 nm light, and an excitation intensity of 12.24 W/cm 2 . The parameter t on was then extracted from resultant simulation data as the time constant of exponentials fit to cytosolic depletion curves and then compared to experimental data. In silico transposition between experiments Apparent t on for the cytosolic depletion method was calculated by fitting an exponential to the cytosolic decay curve. Apparent t on for the colocalization method was calculated as described previously (Glantz et al., 2018); briefly, the correlation coefficient between a GFP-CAAX membrane marker and BcLOV4 was calculated over time for a 4 mm line profile perpendicular to the membrane. The appearance of a membrane virtual marker in the cell was simulated by setting the membrane to a uniform fluorescence intensity and setting the cytoplasm to a uniform intensity equal to 1/10 th that of the membrane. The latter 1/10 th value was calculated from a series of confocal images of GFP-CAAX and represents the cytosolic fluorescence that originates from marker protein that fails to traffic to the membrane. The simulated membrane marker image was convolved with the experimental PSF to account for diffractive effects before downstream analysis. t on could then be calculated by fitting an exponential binding curve to the relationship between correlation coefficient and time.

Virtual cell integration
The system from Equation 1 was input into Virtual Cell and a cell geometry composed of a nucleus, nuclear membrane, cytoplasm, plasma membrane, and extracellular space was also defined. A 3D cell mesh was constructed as described above (see ''Generation of volumetric cell mesh'') and then converted to a stereolithography (,stl) file via custom MATLAB script for import into Virtual Cell. Similarly, the 3D excitation light distribution was computed as described in ''Sensitivity to optical hardware diffraction.'' To describe the volumetric spatiotemporal complexity of laser rastering, the resultant expressions were the scaled summation of multiple offset excitation PSFs, each offset by the spacing between raster points. A ''light'' input element was defined in Virtual Cell and assigned a Cell Reports Methods 2, 100245, July 18, 2022 e12 Article ll OPEN ACCESS value of either 0 when the light was off, or the aforementioned 3D light distribution expression when the light was on. The model could then be run to recover results similar to those of the FEM described herein.
We noted that differences between the Virtual Cell result and the result from our custom system largely stem from the coarser meshing required in Virtual Cell, which necessitates meshing of the extracellular space and is therefore less memory efficient. As a consequence of this coarser meshing, the modeled excitation volume is offset from its experimental location in space. To correct for this offset, the entire cell was shifted along the z axis so that the z-location (relative to the cell) of the simulated excitation volume was re-aligned to the experimental volume.

QUANTIFICATION AND STATISTICAL ANALYSIS
Details of specific statistical analyses for each section, sample sizes, and statistical tests used are given in the STAR Methods and in the corresponding figure legends, but for completeness salient points are also summarized here. Notably, all analysis was done in MATLAB using custom written scripts, which are available with deposited data.

Confocal model accuracy comparison
The performance of the 3D non-linear model created herein were compared to those of a 3D linear model (described in ''Simplified linear model''), a 2D non-linear model, and a 2D linear model. For the 3D non-linear and linear models, mesh generation proceeded as described above. For the 2D non-linear and linear models, the 2D segmented cell frame was directly meshed with the built-in MATLAB mesh generator. The models ran for 100 s for 17 different cells spanning 0.67-10% duty cycles (0.1 s stimulation per 1-15 s period) and were compared by the summed mean-squared error (MSE) between model prediction and experimental result over the entire time-course. Specifically, the MSE refers to error in spatial average cytosolic concentration of BcLOV4 between the model and experimental cell, taken over each time-point in the time-series. If we denote average model-computed cytosolic BcLOV4 concentration at time t as f m (t) and average measured experimental cytosolic BcLOV4 concentration at time t as f e (t), then the MSE represents the summation of (f m (t i ) -f e (t i )) 2 over all timepoints t i in the 100 s experiment (typically every 1 s) divided by the number of timepoints i.
For experimental data, the average over space was calculated by tracking average concentration over time within a 2 mm 3 2 mm square region. This square region was selected based on requirements that it be sufficiently far from the large intracellular voids (i.e., nucleus) and from the membrane edge as to be unaffected by diffractive effects. Based on the PSF of the confocal microscope, the region was >1 mm away from any of these structures. The average fluorescence/concentration timeseries between different intracellular ROIs selected using these criteria tended to be similar (example shown in Figure S5C).
For 3D simulation data, the spatial summation was calculated by first generating 2D slices through the simulation volume at the same z-location as the cell focal plane. For 2D simulation data, this step could be skipped, as resultant simulation images were already 2D slices. The spatial averaging was then done by quantifying average concentration throughout the entire (now 2D) cell cytosol over time. Notably for simulation data, it was not necessary to track fluorescence within a matched ROI, since in the absence of simulated diffractive effects, there were negligible (<1%) subcellular differences in cytosolic BcLOV4 concentration in response to whole-cell stimulation. When analyzing imaging data where subcellular differences in cytosolic depletion may be notable (i.e., for patterned stimulation), comparisons between experimental and model results were done by tracking fluorescence between identical ROIs.
Note that all MSE comparisons were done between 2D images; 3D models were converted to 2D images before MSE analysis. Therefore the metric is equally meaningful for both 2D and 3D models. Statistical comparisons were conducted by paired Wilcoxon signed rank test. Significance was defined as p < 0.05.