Spatially resolved in silico modeling of NKG2D signaling kinetics suggests a key role of NKG2D and Vav1 Co-clustering in generating natural killer cell activation

Natural Killer (NK) cells provide key resistance against viral infections and tumors. A diverse set of activating and inhibitory NK cell receptors (NKRs) interact with cognate ligands presented by target host cells, where integration of dueling signals initiated by the ligand-NKR interactions determines NK cell activation or tolerance. Imaging experiments over decades have shown micron and sub-micron scale spatial clustering of activating and inhibitory NKRs. The mechanistic roles of these clusters in affecting downstream signaling and activation are often unclear. To this end, we developed a predictive in silico framework by combining spatially resolved mechanistic agent based modeling, published TIRF imaging data, and parameter estimation to determine mechanisms by which formation and spatial movements of activating NKG2D microclusters affect early time NKG2D signaling kinetics in a human cell line NKL. We show co-clustering of NKG2D and the guanosine nucleotide exchange factor Vav1 in NKG2D microclusters plays a dominant role over ligand (ULBP3) rebinding in increasing production of phospho-Vav1(pVav1), an activation marker of early NKG2D signaling. The in silico model successfully predicts several scenarios of inhibition of NKG2D signaling and time course of NKG2D spatial clustering over a short (~3 min) interval. Modeling shows the presence of a spatial positive feedback relating formation and centripetal movements of NKG2D microclusters, and pVav1 production offers flexibility towards suppression of activating signals by inhibitory KIR ligands organized in inhomogeneous spatial patterns (e.g., a ring). Our in silico framework marks a major improvement in developing spatiotemporal signaling models with quantitatively estimated model parameters using imaging data.

Introduction Natural Killer (NK) cells are lymphocytes of our innate immune system which provide important immune protection against viral infection and tumors [1,2]. NK cells express a wide range of germ line encoded activating and inhibitory NK cell receptors (NKRs). In humans, activating NKRs include NKG2D and killer Ig-like receptor (KIR) KIR2DS1, and inhibitory NKRs include inhibitory KIRs such as KIR2DL1, and KIR2DL2. NKRs bind to cognate ligands expressed by target cells initiating biochemical, physical, and mechanical modifications within NK cells that culminate into diverse NK cell responses ranging from a neutral response to lysis of target cells to secretion of cytokines [3]. Healthy host cells express polymorphic class I HLA molecules, cognate to a wide range of inhibitory NKRs, and generate tolerant NK cell responses. Whereas, tumor and infected cells downregulate expression of class I HLA molecules or upregulate expressions of ligands cognate to activating NKRs and tip the balance between activating and inhibitory signals toward NK cell activation. NKG2D is one of the best-studied activating NKR. In humans it binds to two families of ligands: one family (MICA and MICB) is akin to MHC class I and the other (ULBP1-3) is related to human proteins that bind to UL16 protein of human cytomegalovirus [4]. Tumor or infected host cells upregulate expressions of NKG2D ligands which contribute to lysis of these cells by NK cell cytotoxic responses [4].
Spatial clustering of NKG2D has been well investigated in confocal [5], total internal reflection fluorescence (TIRF) [6], and super-resolution microscopy experiments [7,8]. Abeyweera et al. [6] using TIRF microscopy experiments reported formation of mobile and immobile microclusters of NKG2D in human NK cell line NKL stimulated by cognate ligand ULBP3 presented on a planar lipid bilayer supported by a glass coverslip. The NKG2D microclusters form at the interface between the NK cell and the lipid bilayer which is also known as the immunological synapse (IS). NKG2D microclusters that are generated at the periphery of IS migrate to the central region of the IS at later times, whereas NKG2D microclusters formed at the central region of the IS remain immobile. The mobility of NKG2D microclusters depends on actin remodeling as treatment by latrunculin, an actin depolymerizing agent, abrogate microcluster movements [6]. Additionally, confocal microscopy experiments show simultaneous localization of NKG2D receptors, Grb2, and Vav1 in the IS in NK cell line NKL upon stimulation by NKG2D ligand MICA [7]. Phosphorylation of Vav1 induces actin remodeling via activation of Rac GTPases [9,10], and thus can regulate motility of NKG2D microclusters.
Spatially resolved computational models have been successfully employed to glean mechanisms that underlie formation, motility, and function of spatial clusters of activating [11] and inhibitory [12,13] NKRs. Kaplan et al. [11] developed a spatially resolved model to investigate hypotheses regarding signal integration of activating NKG2D and inhibitory human KIR2DL/ mouse Ly49 receptors and concluded inhibitory NKRs locally suppress activating signals. Spatial in silico models describing clustering and signaling of inhibitory NKRs elucidated mechanisms giving rise to peptide antagonism for inhibitory KIR2DL2/3 [12]. However, the above in silico models do not quantitatively match or fit microcluster formation and movements of NKRs with that observed in microscopy experiments. As we reason below, in silico modeling of spatial kinetics of NKG2D microclusters is important for gleaning mechanisms and generating improved model predictions. (1) Formation of NKG2D microclusters can increase the production of phosphorylated Vav1 (pVav1) due to increase in the frequency of ligand (e.g., ULPB3) rebinding to NKG2D residing within microclusters, and/or increase in biochemical propensity of signaling reactions when NKG2D molecules are co-clustered with other signaling molecules (e.g., Vav1). (2) Increased production of pVav1 due to clustering of NKG2D can increase centripetal movements of NKG2D microclusters leading to higher spatial aggregation of NKG2D thereby further increasing pVav1 production. This chain of events constitutes a "spatial" positive feedback [14]. Therefore, sizes and spatial separations between NKG2D microclusters could be relevant for affecting downstream signaling.
Motivated by above reasoning and a potential interplay between signaling kinetics and spatial clustering, we developed a framework combining a spatially resolved mechanistic agent based model and published TIRF imaging data assaying spatiotemporal signaling kinetics in NKL cell lines stimulated via NKG2D receptor or NKG2D and KIR2DL2 receptors. The agent based model developed here represents a major improvement over previous modeling efforts in the following aspects: (1) The model is able to quantitatively describe micron scale details of spatial clustering of activating and inhibitory NKRs. (2) A detailed parameter estimation of model parameters is carried out using spatial data. (3) The model includes interplay between spatial clustering and activating NKR signaling kinetics. We investigated three hypotheses (Model 1, Model 2, and Model 3) to show that co-clustering of Vav1 with NKG2D rather than ligand rebinding of ULBP3 is required to produce increased pVav1 production due to the formation of NKG2D microclusters. The presence of the spatial positive feedback allows for an efficient suppression of early time NKG2D signaling by heterogeneously distributed inhibitory HLA-C ligands on target cells.

Model development
We developed a spatially resolved agent based model involving activating NKG2D receptors, inhibitory KIR2DL2 receptors, cognate NKR ligands (ULBP3 and HLA-C), and signaling proteins: Src family kinases (SFKs), Vav1, and phosphatase SHP1. The NKRs and signaling proteins interact via different rules to describe membrane proximal signaling events in human NK cell line NKL. The model includes rules describing biochemical signaling reactions, spatial movements of signaling complexes, and regulation of the spatial movements by biochemical reactions. The rules in the model are designed to quantitatively capture several distinguishing features of spatiotemporal NK cell signaling kinetics observed in previous imaging experiments, namely, a) formation of mobile and immobile NKG2D microclusters upon NKG2D stimulation [6], b) centripetal movements of the mobile NKG2D microclusters from the periphery to the central region of the IS [6,15], c) decrease in the velocity of mobile NKG2D microclusters as those move closer to the central region of the IS [6], and, d) random movements of NKG2D microclusters in the central region of the IS [6].
Our agent based model describes spatiotemporal membrane proximal NK cell signaling in a quasi three-dimensional simulation box representing the interface between NK cell membrane and the plasma membrane of target cell or the lipid bilayer in a TIRF experiment at the IS (Fig 1). The simulation box of area (A = 15μm × 15μm) and thickness z is discretized into chambers of volume l 0 × l 0 ×z where, z = l 0 or z = 2l 0 (l 0 = 0.5μm) for molecules residing in the plasma membrane or the cytosol, respectively. The following spatiotemporal processes occur in the models.
(i) Biochemical signaling reactions. We considered key membrane proximal biochemical signaling reactions involved in NKG2D and KIR2DL2 signaling. The specific roles of signaling proteins (adaptors, kinases, GEFs) considered here in regulating NKG2D signaling have been validated in experiments and modeling over the years. A selection of the pertaining literature is cited below and in Table 1. Similar to Mesecke et al. [16], we take a parsimonious approach in creating the NKG2D signaling network where key signaling reactions are  Table 1. (B) Shows spatial movements considered in the agent based models. The simulation box is divided into chambers of volume l 0 × l 0 × z. A ULBP3 bound NKG2D complex at chamber i moves to its nearest neighbor chamber j with a probability p (m-clus) ji . p (m-clus) ji depends on the number of pVav1 molecules in chamber i and the four nearest neighbor chambers. ULPB3 bound NKG2D complexes in a chamber hop to the nearest neighbor chambers with probabilities p left , p right , p down and p up to implement centripetal movements (see main text). Free protein (receptors, kinases, phosphatases, Vav1/pVav1) molecules hop to next nearest neighbor chambers with probability p diffu implementing diffusive moves. (C) The biochemical signaling induces spatial clustering of NKG2D which in turn increases biochemical signaling-this represents a spatial positive feedback in the models. considered. This approach is widely used for developing "detailed but manageable" signaling kinetic models where a model contains select reactions that pertain to the experimental data the model aims to describe [17]. Below we describe signaling reactions and specific approximations considered in the model. NKG2D receptors bind to cognate ligands (ULBP3) to form NKG2D-ULPB3 complexes. NKG2D homodimers are associated with a pair of DAP10 homodimers in human NK cells. The tyrosine residues in two YINM motifs in a DAP10 homodimer are phosphorylated by Src family kinases [18,19] upon ULBP3 binding. NK cells contain several Src family kinases including Lck, Fyn, Lyn, and Yes [20]. During signaling the DAP10 molecules associated with NKG2D homodimers can be in a variety of partially phosphorylated states where one, two, or three of the total four tyrosine residues are phosphorylated. To reduce the number of agents in the model we approximated partially and fully phosphorylated states of the two DAP10 homodimers by two states, unphosphorylated or fully phosphorylated. In the model, a kinase molecule (SFK) represents multiple Src family kinases that phosphorylate DAP10 upon the formation of the NKG2D-ULBP3 complex. In NK cells, phosphorylated DAP10 becomes available to bind Grb2-Vav1 complex [1,16,21]. In the model, these reactions are approximated by binding of Vav1 to fully phosphorylated DAP10 where Grb2 is not included explicitly. Tyrosine residues in Vav1 have been found to be phosphorylated by SFKs in in vitro assays [22,23]. In NK cells, engagement of adhesion receptor LFA-1 [24] or 2B4 receptors [25] has been reported to lead to phosphorylation of tyrosine residues in Vav1. In an experimentally validated in silico NKG2D signaling model Mesceke et al. [16] considered phosphorylation of tyrosine residues in Vav1 by SFKs. Given the above background we considered DAP10 bound Vav1 is phosphorylated by the SFK in the model. Vav1 phosphorylation is an important event during early time NK cell signaling as pVav1 leads to actin polymerization and degranulation in NK cells [2,26]. Inhibitory KIR2DL2 receptors bind to cognate ligands (HLA-C) and tyrosine residues in immunoreceptor tyrosine based inhibitory motifs (ITIMs) associated with the cytoplasmic part of KIR2DL2 are phosphorylated by the SFKs [27]. SFKs have been implicated in phosphorylation of tyrosine residues on ITIMs [27,28], however, the precise mechanisms are not clear [29]. Phosphorylation of tyrosine residues in ITIMs results in recruitment of SHP-1 [29,30] which lead to dephosphorylation of pVav1 [29,31]. We assume the catalytic domains of SFK phosphorylate the ITIMs associated with KIR2DL2. The model assumes two states of ITIM phosphorylation (unphosphorylated and fully phosphorylated) and the phosphatase SHP-1 binds to fully phosphorylated ITIMs. ITIM bound SHP-1 dephosphorylates pVav1 via enzymatic reactions. Unbinding of ligands from cognate receptors (NKG2D or KIR2DL2) dissociates the signaling complexes completely in the model-this step is included to implement kinetic proofreading [32,33]. In addition, there are first order dephosphorylation reactions for pVav1 representing dephosphorylation of phospho-tyrosines by phosphatases other than SHP1 [34]. Some of the biochemical signaling reactions were not included in the model to keep the model simple and to stay focused on questions of interest, which is further discussed in the Limitations of the model section at the end.
(ii) Spatial movements. We modeled movements of NKG2D-ULBP3 complexes for formation of NKG2D microclusters and centripetal movements of NKG2D microclusters towards the center of the IS. These movements are assumed in the model to be dependent on actin remodeling which is regulated by signaling products such as pVav1 [35]. The above movements are implemented by hops to nearest neighbor chambers occurring with specific probabilities. The velocity of the microclusters decreases as those move closer to the center. In addition, unbound molecules of receptor, ligand, and signaling proteins perform diffusive random movements in the model. Further details are provided in Table 1, Materials and Methods section, and the Supplementary Material. Binding: k (NKG2D-ULBP3) on n NKG2D × n ULPB3 k on estimated by PSO. Range fixed by using K D = 4μM for ULBP3 [36] and k off = 0.023s -1 (for MICA) [36] k (NKG2D-ULBP3) on : k off fixed to 0.023s -1 , also close to the measured k off for ULPB1 [37].
k (DAP10-SFK) phospho : 0.01-10 s -1 Upper limit is chosen to be 1000 × the reported value to match the kinetics at seconds scale in the simulations.

Lck
Set to 698 molecules/(μm) 2 based on the estimated number of NKG2D in a single NKL cell of diameter 18μm [16].

Vav1
Set to 114 molecules/(μm) 3 based on the estimated number of NKG2D in a single NKL cell of diameter 18μm [16].

SHP1
Set to 2090 molecules/(μm) 3 based on the estimated number of NKG2D in a single NKL cell of diameter 18μm [16].

ULBP3
Parameter estimated by PSO. 1075-3468 molecules in the simulation box.
respectively. The rules capture the interplay between biochemical signaling reactions and regulation of specific spatial movements by signaling reactions which in turn affects signaling reactions. For NKG2D signaling the above interplay represents a positive feedback [14]. Hypotheses considered. We considered three hypotheses encoded in models Model 1, Model 2, and Model 3. The models probe functional outcomes of different types of NKG2D clustering and interplay between NKG2D microclusters and biochemical signaling reactions. In Model 1, mobile NKG2D microclusters and NKG2D bound signaling protein molecules move toward the central region of the IS. In Model 2, NKG2D complexes and membrane proximal Vav1 molecules move simultaneously toward the central region of the IS. This Vav1 species in the model could potentially represent Vav1 molecules that are recruited to the plasma membrane via SOS1 [51]. Because of the above rule, there is a strong co-clustering of NKG2D and Vav1 in Model 2. In Model 3, we studied outcomes of the absence of the positive feedback between NKG2D microcluster movement and signaling reactions. In Model 3, the rate of centripetal movements of mobile NKG2D microclusters is independent of pVav1 abundances. Model 3 contains co-clustering of NKG2D and Vav1 as in Model 2. The models are summarized in Table 2.

Model simulation and parameter estimation
The simulations are carried out in a quasi-three-dimensional simulation box representing a thin junction between NK cell and the supported lipid bilayer in TIRF experiments. The simulation box has an area of 15 × 15 μm 2 and a depth of z, and is divided into small cubic chambers of size (l 0 × l 0 × z), where l 0 = 0.5μm and z = l 0 (for molecules residing in membrane) or 2l 0 (for cytosolic molecules). The molecules are well-mixed in each chamber and molecules in a chamber hop to next nearest neighbor chambers with specific rates to produce diffusive or centripetal movements, or movements leading to microcluster formation. The kinetics of the system is simulated using a kinetic Monte Carlo (kMC) approach via a freely available simulator SPPARKS (https://spparks.sandia.gov/). The kMC simulation includes intrinsic noise fluctuations in biochemical reactions as well as in spatial movements. The list of the processes, and their propensities are listed in Table 1. The copy numbers of most of the signaling proteins in the simulation box are estimated using available measured concentrations for NKLs (Table 1). However, values of many of the model parameters in the cellular environment are unknown, and we estimated these parameters by a parameter fitting scheme that reproduced the spatial pattern of NKG2D receptors measured in TIRF experiments [6]. The spatial patterns of NKG2D in our simulation and TIRF imaging data are quantified using mean values, variances, and a two-point correlation function computed from density profiles for NKG2D. Similar variables are widely used in statistical physics [52] to quantify spatial patterns. The Euclidean distance between dimensionless forms of the above variables in TIRF images and the agent based model is used to create a cost function which is minimized by particle swarm optimization (PSO) to estimate model

Multiple models quantitatively describe kinetics of NKG2D microclusters
Spatiotemporal signaling kinetics of NKG2D and associated signaling proteins was simulated in Model 1 and Model 2. NKG2D, ULBP3, SFK, and Vav1 molecules were distributed homogeneously in space in the simulation box at the beginning of the simulation at t = 0. Binding of NKG2D and ULBP3 initiates a series of biochemical signaling reactions (Table 1) leading to phosphorylation of Vav1 in the simulations. The production of pVav1 molecules initiates NKG2D microcluster formation and centripetal movements of the NKG2D microclusters (S1 and S2 Movies). Both models fit spatial distribution of NKG2D in TIRF experiments at t = 1min, in particular, at length scales � 1μm reasonably well (Fig 2). The variances and the two-point correlation functions calculated using spatial distributions of NKG2D in model simulations agreed well with that calculated from the TIRF imaging data (Figs 2 and S1). The estimated best-fit parameters for the models are shown in Table 3. Many parameter values show order of magnitude differences between Model 1 and Model 2, e.g., binding-unbinding rates of SFK to DAP10 or of Vav1 to pDAP10. The reason for this difference can be explained as follows. In Model 2, co-clustering of Vav1 and NKG2D increases propensities of reactions that influence pVav1 production and consequently the formation and motility of NKG2D microclusters. Therefore, different values for reaction rates (e.g., binding/unbinding rate of Vav1 to pDAP10) regulating microcluster kinetics are chosen as optimal parameter values in Model 1 to compensate for the absence of the increase in reaction propensities due to Vav1-NKG2D co-clustering in the model. The computation of uncertainties in the estimated parameter values showed that about 2/3 of the total number of parameters are estimated well, e.g., (standard deviation)/(estimated value) < 3. The procedure for estimation of uncertainty is described in detail in the Materials and Methods section, and potential causes behind poor parameter estimations of 1/3 of the parameters are discussed in the Discussion section. Next, we assessed the utility of the parameter estimation scheme for generating predictions at future time points (Figs 3 and S2) and providing mechanistic insights (next section). The best-fit models were used for predicting spatial patterns of NKG2D at later times t > 1 min. Both models predict spatial distribution of NKG2D in TIRF imaging data until t = 3 min reasonably well (Figs 3 and S2). However, model predictions deviate from the C(r,t) data calculated for the TIRF images at t = 4 min (S3 Fig). This disagreement is likely due to change in the organization of the NKG2D microclusters caused by the spreading of the NK cells on the supported lipid bilayer which is not included in our agent based models.

Co-clustering of NKG2D and Vav1 is required to increase the production of pVav1 due to the formation of NKG2D microclusters
We investigated the mechanistic role of formation and centripetal movement of NKG2D microclusters in increasing pVav1 production. The average lifetime of a NKG2D-ULBP3 complex within NKG2D microclusters can increase because of the increase in the frequency of ULBP3 rebinding due to higher density of NKG2D molecules in microclusters. The increased ULBP3 rebinding could elevate abundances of NKG2D-ULBP3 complexes leading to higher pVav1 production in Model 1 and Model 2. In addition, in Model 2, the increase in reaction propensities due to co-clustering of NKG2D and Vav1 can enhance pVav1 production. In order to evaluate the roles of ULBP3 rebinding and NKG2D and Vav1 co-clustering in increasing pVav1 production we compared kinetics of pVav1 production in Model 1 and Model 2 under two conditions: Case A: NKG2D is not allowed to form microclusters or perform centripetal movements, and, Vav1 does not co-cluster with NKG2D. This case represents NKG2D stimulation in experiments where NKG2D microcluster formation and migration is blocked by application of drugs inhibiting actin polymerization [14]. Case B: Spatial aggregations of NKG2D and Vav1 occur according to the model rules. Simulations for Case A result in spatially homogeneous distribution of NKG2D molecules in both models. Our simulations for Model 1 show that abundances of pVav1 at t = 1 min for a range of ULBP3 doses have negligible differences between Case A and B (Fig 4A). The kinetics of pVav1 production for a particular ULBP3 dose also shows almost no difference between Case A and B (S4A Fig). In contrast, pVav1 abundances decrease in Case A compared to Case B in Model 2 for a range of ULBP3 doses (Figs 4B and S4B). The magnitude of this decrease increases with the increasing ULBP3 dose (Fig 4B). These results suggest that co-clustering of NKG2D and Vav1 in Model 2 could be important for increased pVav1 production in the model. However, ULBP3 rebinding could also help in elevating pVav1 production in Model 2. Therefore, we further quantified the contribution ULBP3 rebinding in increasing the average number of NKG2D-ULBP3 complexes in Model 2. We followed an approach in ref. [53] for this quantification, wherein decay kinetics of an initial fixed number of receptor-ligand complex is studied in the presence of immobile spatially distributed receptors. The simulations start with a fraction of receptors bound to ligands and no free ligands; free ligands created by dissociation of the receptor-ligand complex at a rate k off can diffuse and rebind to the receptors. In the absence of any ligand rebinding, the number of receptor-ligand complex decays exponentially as / exp(-k off t). The presence of rebinding produces a slower and a non-exponential decay of the number of receptor ligand complex with higher numbers of receptorligand complex remaining in the system at longer times (t� 1/k off ) compared to an exponential decay as exp(-k off t). We evaluated the decay kinetics of an initial number of NKG2D-ULBP3 complex where NKG2D-ULBP3 binding-unbinding reactions were the only reactions present in simulations. The NKG2D molecules, held fixed in space, were distributed uniformly and randomly or in a spatially clustered configuration obtained from our simulations for Model 1 or Model 2 at t = 1 min (Fig 2C). The decay kinetics shows non-exponential decay for both the cases and a small increase in the number of NKG2D-ULBP3 complex ( Fig  4C and 4D) when NKG2D are clustered. The above results reveal that the effect of ULBP3 rebinding in our models is almost equally strong when NKG2D are distributed uniformly randomly or in microclusters. Thus, increased ULBP3 rebinding in NKG2D microclusters does not play a substantial role in increasing pVav1 production.
The increase in the production of pVav1 induced by NKG2D microclusters in Model 2 ( Fig 4B) is consistent with experiments by Endt et al. [14], where blocking cytoskeletal movements by an actin polymerization inhibiting drug cytochalasin D produces a substantial decrease in pVav1 in NKLs stimulated by NKG2D ligand MICA. Actin polymerization induces microcluster formation and centripetal movements of NKG2D as treatment with actin depolymerization agent latrunculin abrogated NKG2D microcluster formation and movements in NKLs [6]. Next, we used Model 2 for deciphering mechanisms of signal integration.

pVav1 dependent centripetal movements of NKG2D are abrogated by inhibitory KIR2DL2 signaling
We investigated inhibition of NKG2D signaling by inhibitory KIR2DL2 receptors in Model 2. The KIR2DL2 receptors were distributed in the simulation box following TIRF imaging data in ref. [6] (S5A Fig). The two-point correlation function calculated from spatial distribution of KIR2DL2 at t = 0 in our simulations agrees well with that of the TIRF imaging data (S5C Fig). The TIRF experiments show that KIR2DL2 microclusters are present in higher numbers at the periphery compared to the central region of the IS. In addition, the spatial organization of these microclusters does not change appreciably between pre-and~30s post-stimulation by HLA-C ligands (S5B Fig). The changes in KIR2DL2 microcluster distributions beyond 30s occurred presumably due to rapid retraction of NK cells from the supported lipid bilayer. Since we did not model NK cell retraction in our models, we assumed KIR2DL2 microclusters to be stationary in our simulations. The chambers where the number of KIR2DL2 exceeded a threshold value were considered to be parts of KIR2DL2 microclusters in simulations. In our models, KIR2DL2 microclusters colocalize with SFK which is supported by previous experiments [54]. NKG2D, ULPB3, HLA-C, SFK, SHP1, and Vav1 were distributed homogeneously in the simulation box at t = 0. Some of the SFK molecules were co-clustered with KIR2DL2 which remained immobile for the duration of the simulations (details in Materials and Methods section). In the simulations, SHP1 recruited by the pITIMs in KIR2DL2-HLA-C complexes dephosphorylate pVav1 (bound to NKG2D complex or free) residing in the same spatial location or the same chamber. Thus, the production of pVav1 reduces substantially in the presence of inhibitory KIR2DL2 signaling (Fig 5B). The decrease in pVav1 abundance also slows down centripetal movements of the NKG2D microclusters resulting in a more spread-out spatial distribution of NKG2D in simulations (Fig 5A and S3 Movie). We quantified the reduction of centripetal migration of NKG2D by calculating the increase in the number of NKG2D molecules in a region enclosing the center of the IS (Fig 5C) over a time period. This decrease in the centripetal movement of NKG2D is qualitatively in agreement with experiments by Abeyweera et al. [6].

Interplay between Vav1 phosophorylation and centripetal movements is required for efficient inhibition of NKG2D signaling by inhibitory KIR2DL2
We employed Model 2 and Model 3 to study the role of the interplay between centripetal movements of NKG2D microclusters and pVav1 in early time integration of NKG2D and KIR2DL2 signals. In Model 3, centripetal movements of NKG2D and Vav1 are independent of pVav1 abundances, thus, the decrease in pVav1 abundances due to inhibitory signaling does not affect accumulation of NKG2D at the central region of the IS (S4 Movie and Fig 6B). We carried out simulations for two scenarios where inhibitory HLA-C ligands were distributed in the simulation box homogeneously (Case I) or in a ring pattern (Case II) (Fig 6A) devoid of HLA-C molecules at the center of the ring. The ring pattern could represent HLA molecules on engineered surfaces [55] or on target cells. Almeida et al. [54] found HLA-C molecules to be organized in a ring pattern on target cells after NK cells are incubated with target cells for 10mins. The pattern formation in part is initiated by binding of inhibitory KIR2DL2 and HLA-C, and KIR2DL2 co-localizes with HLA-C in such clusters. Simulations were carried out with HLA-C and ULBP3 in Model 2 and Model 3. When HLA-C are distributed homogeneously (Case I), the total number of pVav1 at t = 1 min post NKG2D stimulation decreases substantially in both models as the number of HLA-C increased about 4 -fold (1000 to 4000 in the simulation box) (Fig 6D). The pVav1 abundances are slightly higher in Model 3 than Model 2 owing to increased NKG2D clustering in the former model. This decrease in pVav1 with increasing HLA-C in the models is qualitatively similar to the large decrease in the percentage of lysis of human NK cell line YTS in cytotoxicity assays reported in ref. [54] as the abundance of HLA-C on target cells increased about 4 -fold. Next, we investigated variation of pVav1 with increasing HLA-C dose in Model 2 and Model 3 when HLA-C molecules were distributed in a ring pattern (Case II). The simulations show smaller decrease in the number of pVav1 in Model 3 compared to Model 2 as the number of HLA-C is increased (Fig 6C). This result can be explained as follows. In Model 3, NKG2D microclusters, regardless of the pVav1 concentrations, accumulate at the center of the IS devoid of HLA-C and produce pVav1 even at high HLA-C concentrations, whereas, in Model 2 the lower amount of pVav1 at high HLA-C concentrations decreases centripetal movements of NKG2D and prevents accumulation of NKG2D at the center of the IS (S6 Fig). We also carried out the above comparison when KIR2DL2 molecules are co-localized with HLA-C in the ring pattern to represent HLA-C clustering in target cells interacting with NK cells [54]. The results (S7 Fig) were similar to that in Fig 6C. Therefore, the HLA-C ligands are able to suppress the production of pVav1 more efficiently in Model 2 compared to Model 3.

Discussion
Spatial organization of activating human NKRs such as CD16 [56], NKp46 [57], KIR2DS1 [58], and NKG2D [5][6][7], and associated signaling proteins has been observed in NK cells stimulated by cognate activating ligands. Microcluster formation by NKRs can increase the apparent lifetime of NKR-ligand complexes because of higher frequency of ligand rebinding within a microcluster. In addition, co-clustering of NKRs with signaling proteins in NKR microclusters can increase biochemical propensities of signaling reactions due to elevated local concentrations of reacting proteins. Either of the above mechanisms could increase production of activated (e.g., tyrosine phosphorylated) signaling proteins. However, when the contribution due to ligand rebinding is not substantial, co-clustering of NKR and other signaling proteins plays a dominant role in increasing downstream signaling. We developed a predictive in silico framework to quantify relative roles of the above mechanisms for early (~few mins) NKG2D signaling kinetics to demonstrate that co-clustering of Vav1 with NKG2D plays a more dominant role over ULBP3 rebinding in increasing activating signals. Since NKG2D binds with ULPB3 with a half-life (ln(2)/k off ) of~30 seconds and the average number of homogeneously distributed NKG2D molecules in the area occupied by a typical NKG2D microcluster is about two molecules, the higher numbers of NKG2D in microclusters do not result in a large increase in abundances of NKG2D-ULBP3 complexes due to ligand rebinding. Our simulations confirm the above explanation. Co-clustering of signaling proteins could play an important role for increasing downstream signaling for other activating NKRs as well. The density of CD16 on primary NK cells is about 200 times larger than that of NKG2D [59] and the lifetime of a CD16-ligand complex is similar (e.g., k off <0.01 s -1 for IgG [60]) or larger than that of NKG2D-ULBP3, therefore, the effect of co-clustering of signaling proteins with CD16 could be more relevant for increasing downstream signaling than ligand rebinding. Our NKG2D signaling model could describe spatiotemporal clustering by other activating NKRs such as CD16 qualitatively, as the formation and movements of another activating NKR microclusters could arise from similar interplay between pVav1 kinetics and cytoskeletal reorganizations modeled here, however, the time scales of formation and centripetal movements of microclusters of activating NKRs would depend on the time scales of pVav1 production. Since, the signaling reactions leading to pVav1 production can be different between NKG2D and another activating NKR (e.g., CD16, KIR2DS1), the spatiotemporal kinetics of NKG2D and other activating NKRs are likely be quantitatively different.
Our agent based model (Model 2) is successful in generating predictions regarding inhibition of NKG2D signaling by inhibitory KIR2DL2 signaling or by drugs inhibiting actin polymerization. The model qualitatively captures the slowing down of the NKG2D microcluster movements in the presence of inhibitory KIR2DL2 signaling (Fig 5A and 5C) and the decrease in pVav1 (Fig 5B) in response to an actin polymerization inhibiting drug cytochalasin D in NKLs stimulated by NKG2D ligand MICA [14]. The model also quantitatively predicts spatial pattern of NKG2D in TIRF imaging at later time points, not included in model training, reasonably well. The model also reasonably predicted (S8 Fig) early (� 70 secs) clustering of NKG2D in NKLs stimulated by MICA on target Daudi cells in confocal experiments reported by Brown et al. [7]. This provides confidence on the applicability of our model and parameter estimation for describing other systems. The model can be further expanded to investigate mechanisms of interplay between pairs of activating and inhibitory NKRs stimulated by activating ligands and different organizations of HLA-C ligands with peptides. Adaptable spatial patterns of antigens presented by DNA origami structures have been recently used to manipulate responses in B-cells [61] and CAR T-cells [62]. The spatial model developed here can serve as a screening tool for identifying spatial patterns of antigens optimal for generating specific lymphocyte responses.
In order to investigate the utility of our model to predict further downstream NK cell responses, we applied our model to predict lysis of target cells treated with DNA polymerase inhibitor aphidicolin in NK cell cytotoxic assays by Gasser et al. [63]. Gasser et al. [63] observed aphidicolin treatment moderately increases the expression of NKG2D ligands in mouse T-cell blast target cells and these target cells showed increased lysis compared to their untreated counterparts. An extension of our signaling model qualitatively predicts (S3 Text and S16 Fig) the differences in lysis of aphidicolin treated and untreated target cells for a range of effector cell:target cell ratios as reported by Gasser et al. [63]. This application could point to a potential way to combine early time signaling model with population kinetic models describing later time NK cell responses.
We studied the role of the interplay between creation of NKG2D microclusters, centripetal movement of the microclusters, and early time signaling using an alternate model (Model 3) where centripetal movement of NKG2D microclusters is independent of pVav1. Our simulations show that pVav1 dependent centripetal movements allow for a more efficient suppression of pVav1 in the presence of inhibitory KIR signaling, in particular, when inhibitory ligands (HLA-C) are clustered in a ring formation in the target cell membrane. HLA-C molecules have been found to organize in ring or multifocal patterns on target cells [54], the presence of the spatial feedback in NKG2D microclusters and pVav1 creates an efficient and flexible maneuver to inhibit signaling when inhibitory ligands are present in spatially inhomogeneous patterns.
A unique aspect of our modeling framework is estimation of model parameters that best fit spatial patterns of NKG2D measured in TIRF imaging experiments. Many parameters in our model were not measured in experiments previously or are difficult to measure because they described coarse-grained processes (e.g., parameters in describing centripetal force). In addition, reported parameter values measured under specific conditions (e.g., outside the cell) might not reflect its magnitude within the cell. Therefore, we carried out an estimation of model parameters using parallel codes that produced quantitative agreement between spatial patterns of NKG2D in simulation and imaging experiments. We used two-point correlation functions to quantify spatial organization of NKG2D microclusters which is widely used in statistical physics and materials science to quantify spatial patterns. The parameter estimation step provides the following benefits: (i) it allows models to capture unique aspects of kinetics of NKG2D microclusters, e.g., centripetal movement with decreasing mobility at the center of the IS. (ii) Provides parameter ranges where the roles of ligand rebinding and co-clustering of NKG2D and Vav1 can be compared. Since both the processes depend on parameter values, these should be compared in parameter ranges that are able to reproduce experimental data. Receptor clustering in the IS has been widely investigated in computational models in T- [64][65][66][67][68][69][70], B- [71], and NK-cells [11,12], however, most of these models described the spatial patterns qualitatively. Some of the modeling quantitatively described receptor patterns in NK- [72] and T- [73] cells, however, did not couple receptor clustering with signaling kinetics. A major difficulty in parameter estimation in spatial models is computationally intensive nature of such calculations, which we address here using parallel computation. However, as we discuss below there are several areas where this framework requires further improvement.

Limitations of the work
The model did not include several processes or made approximations to keep the model simple and focused on the questions of interest addressed here. We did not include integrin receptor LFA1 and its ligand ICAM-1 in our models for simplicity and focused our study to NKG2D and KIR2DL1 signaling and clustering. Larger size LFA1-ICAM-1 complexes (~40nm) segregate from the smaller sized (~10-15nm) NKG2D-MICA or KIR2DL1 molecules due to size based exclusion interactions [74]. Usually, LFA1-ICAM1 molecules reside in the periphery of the IS surrounding NKG2D [74,75]. However, integrin receptor signaling contributes towards NK cell signaling [76][77][78] and spatial clustering of LFA-ICAM1 will likely contribute towards early time NK signaling and its regulation of NKG2D spatiotemporal signaling kinetics. Similarly, proximal signaling proteins such as PI3K [79] were not included in the model to keep the model simple and focused on the questions of interest addressed here. The models also do not include NKG2D degradation [80] which tends to occur later in the signaling (>15 mins) [81] and could underlie the lower amount of phosphorylated ITAM phosphorylation in the central region in the TIRF experiments in Abeyweera et al. [6]. We represented multiple SFKs (e.g., Lck, Fyn, Lyn, Yes) in NK cells by a single species for simplifying models. There are differences in substrate specificities [82] for SFKs in T cells. It will be straightforward to extend our models to include multiple SFKs if similar differences are found in NK cells. Furthermore, our model did not include spreading of the NK cell surface on the glass slide. Upon NKG2D stimulation the NK cell spreads at later times (� 4min). Due to this reason our model can describe early time (0-3 min) NKG2D spatiotemporal signaling kinetics where NK cell spreading is negligible but is unable to quantitatively predict changes in NKG2D spatial patterns due to spreading of NK cell surface on the glass slide at times � 4min (S3 Fig) or late time signaling events such as secretion of lytic granules. The model also does not include retraction of the NK cell surface from the bilayer within minutes (~5 mins) induced by stimulation of inhibitory KIR2DL2 in the TIRF experiments [6]. The retraction helps to self-limit the formation of NK cell-target cell conjugates when inhibitory signals dominate [5,6].
One the computation side, several parameters showed large error bounds implying these parameters can be changed substantially but will make small or no changes to the mean cluster sizes and two-point correlation functions. A potential solution to address this challenge could be introducing weight factors for combining mean values and two-point correlation functions (or second moments) following a systematic framework such as generalized method of moments [83,84].

Spatial kinetic Monte Carlo simulation
We used the software package SPPARKS (https://spparks.sandia.gov/) to simulate reactions and particle hopping moves described below. We chose Gibson-Bruck implementation of Gillespie algorithm [85] within SPPARKS to perform the kinetic Monte Carlo simulation.

Biochemical reactions
Molecules in each chamber of volume (l 0 × l 0 × z) where l 0 = 0.5μm and z = l 0 (for plasma membrane bound molecules) or z = 2l 0 (for cytosolic molecules) are taken to be well-mixed. Stochastic biochemical reactions in individual chambers are simulated using reaction propensities listed in Table 1 [86]. Therefore, the volume factor (v recep-lig ) used for converting the unit of binding rate (k on ) from (μM) -1 s -1 to propensity (/ k on /v recep-lig ) in s -1 in a chamber is taken as, v recep-lig = l 0 × l 0 × d recep-lig . Similarly, the volume factor v recep-SFK used for binding of NK cell receptors with SFKs is set to, v recep-SFK = l 0 × l 0 × d recep-SFK , where d recep-SFK~1 0 nm [16]. The volume factor v cytosol for cytosolic molecules binding with NK cell plasma membrane bound complexes is taken as, v cytosol = (l 0 ×l 0 ×2l 0 ). In Model 2, the volume factor for binding of plasma membrane associated Vav1 molecules with other plasma membrane residing molecules is taken as, v Vav1-plasma = (l 0 ×l 0 ×l 0 ).

Microcluster formation
A pVav1 dependent processes is implemented to generate microclusters of NKG2D. A "potential" function E i is associated with any chamber i, is given by, where, x i is the number of pVav1 molecules in the i th chamber and {j1, j2, j3, j4} denote its four nearest neighbors. The probability of moving NKG2D complexes from the i th chamber to its nearest neighbor j is given by, p mÀ clus ; thus, when a nearest neighbor j has a potential lower than the i th chamber, i.e., E j < E i , NKG2D molecules move to the j th chamber. We estimate β using PSO.

Diffusive movements
Unbound molecules of NKRs, their cognate ligands, SFK, SHP1, and Vav1 move diffusively. The diffusive movements are simulated by hopping moves of molecules to the nearest neighbor chambers with probability p diffu / D/(l 0 ) 2 . The propensities for these movements are given in Table 1. Periodic boundary conditions in the x-y plane are chosen for freely diffusing molecules. The propensity for diffusive movements of cytosolic molecules is about >1000 fold larger than that of the plasma membrane bound molecules and that of most reactions. This is because of the larger value of the diffusion constant (~10 μm 2 /s) [87] and the larger average number of cytosolic molecules (~30-200 molecules) in a chamber. Therefore, a substantial proportion of the Monte Carlo moves are spent on diffusive moves of cytosolic molecules which make the run time of the simulation long (wall time~10 hrs on a single 3.0 GHz AMD EPYC processor). We approximated these fast diffusive moves by not executing explicit diffusion moves for the cytosolic molecules but by homogenizing cytosolic molecules in the simulation box at time intervals proportional to the average time a cytosolic molecule will need to diffuse the length of the simulation box. We checked the validity of this approximation by comparing simulations incorporating the above approximation with those containing explicit diffusive movements for cytosolic molecules. We did not find appreciable differences between the two (S9 Fig).

Microcluster movements
Any ULBP3 bound NKG2D receptor complex is assumed to be a part of a microcluster. All molecules of a specific ULBP3-NKG2D complex, e.g., ULBP3-NKG2D-SFK, in a chamber are moved to one of the four nearest neighboring chambers in a single microcluster hopping move. The probabilities of the hops to the neighboring chambers are given by p ! hop � (p right , p left , p up , p down ) and the corresponding propensities are calculated by multiplying the probabilities by a rate k (cluster-move) . p ! hop at a position (x,y) in the x-y plane is given by, , where, the four corners and the center of the simulation box in the x-y plane are at (0, 0), (0, 2R), (2R,0), (2R, 2R), and (R,R), respectively. The variable s depends on the total number of pVav1 (or [pVav1] T ) in the simulation box and is given by, The hopping probability ! p hop is composed of two parts: one generates centripetal movements [88,89] and the other produces random Brownian movements. The weight factor 0 � w� 1 determines the relative proportions of centripetal and random components in p ! hop , i.e., the smaller the w, the stronger is the bias toward random movements. The velocities for centripetal movements along x and y directions are given by, Therefore, the magnitude of the centripetal velocity v r is given by, v r ¼ sw ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi : v r decreases monotonically from the periphery to the center of the simulation box. v r = 0 at the center (R,R) and is maximum ( ffi ffi ffi 2 p sw) at the corners of the simulation box. In the absence of any pVav1 in the simulation box, the centripetal movements cease to exist, i.e., v r = 0. Velocity of F-actin retrograde flow decreases monotonically from the periphery to the center in YTS NK cell lines stimulated by surface coated anti-NKG2D [15]. We assume that NKG2D microcluster movements are guided by F-actin, similar to that to TCR microclusters in antigen stimulated Jurkat T cells [90].
The rate, k (cluster-move) determines the magnitude of the centripetal velocity and the diffusion constant of the Brownian movements of microclusters. A movie of movements of a single microcluster is shown in the supplementary material (S5 Movie).

Excluded volume interactions
We impose an upper bound (N thres ) to the number of molecules that can reside in a chamber to incorporate excluded volume interactions. Molecules are allowed to hop to a chamber if the number of molecules does not increase beyond N thres due to the move. A separate upper bound for the number of NKG2D in a chamber is implemented to prevent aggregation of NKG2D to a small number of very high density NKG2D microclusters in the central region.

Influx of NKG2D
An influx of NKG2D molecules in the simulation box from the boundary with a rate k (influx) is implemented by adding an outer rim of thickness 2l 0 at the boundary of the simulation box. The chambers in the rim are occupied by NKG2D and ULBP3 which undergo NKG2D-ULPB3 binding-unbinding reactions. The NKG2D-ULBP3 complexes do not participate in further downstream reactions inside the outer rim. The bound NKG2D-ULPB3 complexes enter the simulation box by hopping moves that occur at rate k (influx) .

Initial configuration
Unbound molecules of NKG2D, SFK, Vav1, KIR2DL2, and SHP1 were distributed homogeneously in the simulation box where any chamber contained the same number of molecules for a particular protein. The numbers of molecules in each chamber for the above protein species are calculated using the values shown in Table 1 assuming proteins are distributed homogeneously in the plasma membrane or the cytosol of an NK cell. The number of ULBP3 in the simulation box is estimated via PSO, where, a chamber at t = 0 is populated with 3 molecules (or left unpopulated) of ULBP3 with a probability f (or 1-f). f is estimated in PSO.

Particle swarm optimization
We performed asynchronous particle swarm optimization (PSO) algorithm with constraints to optimize a cost function to an estimate model parameters. The parameters were varied as powers of 10 in PSO. The parameters used in the PSO are the following: particle velocity scaling factor, ω = 0.5, coefficient weighting a particles best-known position, φ p = 2.5, and coefficient weighting the swarm's best-known position, φ g = 1.5. A swarm size of 200 particles is considered. The maximum iteration limit is assigned to 100. The computation takes about 48h on 300 parallel 3.0 GHz AMD EPYC CPUs. A similar set of PSO parameters were used to estimate parameters in a spatial model describing formation of bacterial biofilms in three dimensions [91]. The construction of the cost function is described below.

Construction of the cost function
We extracted intensity profile I({x i ,y i },t) of NKG2D-DAP10 using TIRF images (details in Extraction of image data), where {x i ,y i } denote the centers of the chambers in the x-y plane in the simulation box. Our parameter estimation scheme minimizes a cost function that measures the Euclidean distance between variables quantifying statistical properties of spatial patterns of NKG2D in TIRF experiments and simulations from our agent based models. We computed mean, variance, and two-point correlation function [52,91], The summation over r x and r y indicates average of all neighbors of (x i ,y i ) separated by a distance r, i.e., r x 2 + r y 2 = r 2 . We used periodic boundary conditions for the calculation of C S (r,t) for r � L/2, where L is the length of the simulation box. For the TIRF imaging data, S({x i ,y i },t) = I(x i ,y i ,t)/I max (t), and for our model simulations, S({x i ,y i }, t) = n({x i ,y i }, t; θ)/n max (t), where n (x i ,y i ,t; θ) denotes the number of NKG2D molecules in a chamber centered at (x i ,y i ) in the x-y plane. θ represents model parameters listed in Table 1. I max (t) and n max (t) denote the maximum values of image intensity and NKG2D number in the chambers in the imaging data and in the simulation, respectively. The quantities in Eq (2) depend on the initial state of the system at t = 0 in simulations or experiments and on intrinsic noise fluctuations that arise during time evolution of the system (S10 Fig). Averages of these quantities (denoted by � � � h i) over ensembles of configurations at time t arising from different initial states and stochastic trajectories with intrinsic noise fluctuations are traditionally used in physics [52] to characterize spatial patterns. However, it is challenging to perform such ensemble averages in our situation because, (i) few replicates of experimental data are usually available, and (ii) computational cost of performing ensemble averaging within PSO. In order to circumvent these issues, we minimized a cost function as a function of the quantities in Eq (2) where the best-fit values of the model parameters depend on the initial state as well as on the intrinsic noise fluctuations. However, we account for effect of the intrinsic noise fluctuations in the parameter estimation as described in the error estimation section.
The cost function is given by, η represents random variables associated with the intrinsic noise fluctuations and n 0 denotes the initial state at t = 0 in the model. The minimization of Eq (3) in our PSO yields a θ = θ min associated with a specific set of random variables η pso-optim and a fixed initial state n 0 , i.e., the minimum value of C cost , C cost | min = C cost (θ min ; η = η pso_otim ; n 0 ).

Quantification of uncertainties in PSO estimation of parameters
We generated configurations {n({x i ,y i },t; θ = θ min )} for the best-fit parameter value θ min for an initial state n 0 , and different sets of intrinsic noise fluctuations (or different η) in our simulations. We computed C cost (θ min ; η, n 0 ) given by Eq (3) for each n({x i ,y i },t; θ min ) and computed the variance, σ C , for the values of C cost evaluated for the ensemble,{n({x i ,y i },t; θ)} (S11 Fig). We reasoned that the range 0� C cost � C cost | min +2σ C , will be scanned by intrinsic noise fluctuations in the model for θ min ; therefore, a configuration n({x i ,y i }, t; θ) where θ 6 ¼ θ min that generates a cost function C cost within above range cannot be separated well from {n({x i ,y i },t; θ min )}. This produces an uncertainty in our estimated parameters θ min . We characterized the uncertainty in θ min by collecting the parameters {θ} that produce cost functions in range 0 to C cost | min +2σ C (S11 Fig). Next, we evaluated if these parameters reside in a single or multiple clusters in the manifold spanned by the parameters-the existence of multiple clusters will indicate the presence of multiple local minima. We computed the number of cluster following a method in Ref. [92] and found a single cluster (S12 Fig). Each parameter was scaled as (θ-θ min )/(θ max − θ min ), where θ min and θ max are the minimum and the maximum parameter values, respectively, to make it dimensionless and to lie between 0 and 1 before applying the algorithm in Ref. [92]. The standard deviations computed for the points in the cluster gives an estimate of the uncertainty in θ min .

Extraction of image data
Intensities for fluorescently labeled molecules in TIRF experiments were extracted from images published in in Ref. 6. The python code used for the extraction is uploaded on github at https://github.com/jayajitdas/NK_signaling_spatial_model. The extracted intensities for fluorescently tagged NKG2D-DAP10 or KIR2DL2 molecules from regions of interest were averaged over multiple pixels to create spatial data with the minimum resolution (~0.5 μm) in our simulation. Further details are provided in S13 Fig . We extracted two-dimensional fluorescence intensity of NKG2D-GFP in single NKL published in Ref. [7] following an image extraction method described in the Materials and Methods section. The NKLs in Ref. [7] were stimulated by MICA ligands on target Daudi/MICA cells and NKG2D-GFP molecules were imaged using confocal microscopy. We chose a region of interest in the extracted image and coarse-grained it to match the minimum length scale (0.5 μm) of our simulation. We used the confocal data at t = 0 (A; top) to create the initial configuration of NKG2D receptors in the model (Model 2). The rest of the parameters are initialized as described in the Materials and Methods sections. We simulated the initial configuration in our model (Model 2) using the best fit parameters (Table 3) we obtained for Model 2 using the TIRF imaging data in Ref. [6].
The spatial organization of NKG2D in the simulation (top panels in (A)-(C)) and experiments (bottom panels in (A)-(C)) at t = 30s and t = 70s were compared using the two-point correla- μ and σ denotes the mean and standard deviation for C cost (θ min ; η, n 0 ) for the corresponding model. C opt represents the optimum (minimum) cost function C cost (θ min ; η = η pso_otim ; n 0 ) estimated by PSO for each model. The uncertainty in our estimated parameters θ min is estimated to lie within the interval 0� C cost � C thres, where C thres is obtained as C thres = C opt +2σ.
(TIF) S12 Fig. Decision Graph to compute clusters according to Density Peak Clustering algorithm (Ref. [92]). The points having relatively large local density (ρ) and high value of δ are considered to be the cluster centers defined in Ref. [92]. For both, Model 1 (A) and Model 2 (B), we observed only one such cluster center to exist implying the presence of a single minimum cost function in the parameter range explored by the PSO. The decision graph was computed for parameters θ in the solution space having cost function, C cost such that 0� C cost � C thres, , where C thres = 0.29 and 0.35 for Model 1 and Model 2, respectively. The parameter d c in the Density Peak Clustering algorithm in Ref. [92] is calculated such that the average number of neighbours is 2% of the total number of points in the data set.  Fig 3A (top). Each value in the .vtk data file corresponds to the intensity of NKG2D-DAP10-mCherry, extracted from TIRF experiments in Ref. [6] at t = 2min, post stimulation by ULPB3, at a x-y grid location in the simulation box. "DIMENSIONS" and "POINT DATA" in the header of the data file represent the dimensions of the simulation box and total number of chambers in the simulation box, respectively. The image is generated using Para-View-5.7.0 which takes .vtk files as inputs.
(VTK) S6 Data. Data for Fig 3A (middle). Each value in the .vtk data file corresponds to the number of NKG2D in each chamber of the simulation box for our Model 1 at t = 2 min, post NKG2D stimulation. "DIMENSIONS" and "POINT DATA" in the header of data file represent the dimensions of the simulation box and total number of chambers in the simulation box, respectively. The image is generated using ParaView-5.7.0 which takes .vtk files as inputs.
(VTK) S7 Data. Data for Fig 3A (bottom). Each value in the .vtk data file corresponds to the number of NKG2D in each chamber of the simulation box for our Model 2 at t = 2 min, post NKG2D stimulation. "DIMENSIONS" and "POINT DATA" in the header of data file represent the dimensions of the simulation box and total number of chambers in the simulation box, respectively. The image is generated using ParaView-5.7.0 which takes .vtk files as inputs.
(VTK) S8 Data. Data for Fig 3B (top). Each value in the .vtk data file corresponds to the intensity of NKG2D-DAP10-mCherry, extracted from TIRF experiments in Ref. [6] at t = 3min, post stimulation by ULPB3, at a x-y grid location in the simulation box. "DIMENSIONS" and "POINT DATA" in the header of the data file represent the dimensions of the simulation box and total number of chambers in the simulation box, respectively. The image is generated using Para-View-5.7.0 which takes .vtk files as inputs.
(VTK) S9 Data. Data for Fig 3B (middle). Each value in the .vtk data file corresponds to the number of NKG2D in each chamber of the simulation box for our Model 1 at t = 3 min, post NKG2D stimulation. "DIMENSIONS" and "POINT DATA" in the header of data file represent the dimensions of the simulation box and total number of chambers in the simulation box, respectively. The image is generated using ParaView-5.7.0 which takes .vtk files as inputs.
(VTK) S10 Data. Data for Fig 3B (bottom). Each value in the .vtk data file gives the number of NKG2D in each chamber of the simulation box for our Model 2 at t = 3 min, post NKG2D stimulation. "DIMENSIONS" and "POINT DATA" in the header of data file represent the dimensions of the simulation box and total number of chambers in the simulation box, respectively. The image is generated using ParaView-5.7.0 which takes .vtk files as inputs.
(VTK) S11 Data. Data for Fig 5A. Each value in the .vtk data file gives the number of NKG2D in each chamber of the simulation box for our Model 2 in the presence of both activating ULBP3 and inhibitory ligands at t = 1 min. "DIMENSIONS" and "POINT DATA" in the header of data file represent the dimensions of the simulation box and total number of chambers in the simulation box, respectively. The image is generated using ParaView-5.7.0 which takes .vtk files as inputs.
(VTK) S12 Data. Data for Fig 6B. Each value in the .vtk data file gives the number of NKG2D in each chamber of the simulation box for our Model 3 in the presence of both activating ULBP3 and inhibitory ligands at t = 1 min. "DIMENSIONS" and "POINT DATA" in the header of data file represent the dimensions of the simulation box and total number of chambers in the simulation box, respectively. The image is generated using ParaView-5.7.0 which takes .vtk files as inputs. (VTK)