An optimal control approach for enhancing natural killer cells' secretion of cytolytic molecules

Natural killer (NK) cells are immune effector cells that can detect and lyse cancer cells. However, NK cell exhaustion, a phenotype characterized by reduced secretion of cytolytic models upon serial stimulation, limits the NK cell's ability to lyse cells. In this work, we investigated in silico strategies that counteract the NK cell's reduced secretion of cytolytic molecules. To accomplish this goal, we constructed a mathematical model that describes the dynamics of the cytolytic molecules granzyme B (GZMB) and perforin-1 (PRF1) and calibrated the model predictions to published experimental data using a Bayesian parameter estimation approach. We applied an information-theoretic approach to perform a global sensitivity analysis, from which we found that the suppression of phosphatase activity maximizes the secretion of GZMB and PRF1. However, simply reducing the phosphatase activity is shown to deplete the cell's intracellular pools of GZMB and PRF1. Thus, we added a synthetic Notch (synNotch) signaling circuit to our baseline model as a method for controlling the secretion of GZMB and PRF1 by inhibiting phosphatase activity and increasing production of GZMB and PRF1. We found that the optimal synNotch system depends on the frequency of NK cell stimulation. For only a few rounds of stimulation, the model predicts that inhibition of phosphatase activity leads to more secreted GZMB and PRF1; however, for many rounds of stimulation, the model reveals that increasing production of the cytolytic molecules is the optimal strategy. In total, we developed a mathematical framework that provides actionable insight into engineering robust NK cells for clinical applications.


INTRODUCTION
Natural killer (NK) cells are innate immune effector cells that protect the host from diseased cells such as virally infected cells and cancer cells. 1,2 In particular, when NK cells engage with these target cells, NK cell stimulatory receptors become activated and mediate killing of the diseased cells. One mechanism for target cell killing is through the secretion of the cytolytic molecules granzyme B (GZMB) and perforin-1 (PRF1). [3][4][5][6] Secretion of these factors is termed "degranulation." Specifically, PRF1 mediates the formation of pores on the target cell membrane, enabling GZMB to infiltrate and induce apoptosis. Although the secretion of cytolytic molecules is mediated by multiple NK cell receptor signaling pathways, 7 including -but not limited to -the natural cytotoxicity receptors (e.g., NKp46), 2B4 (CD244), and DNAM-1 (CD226), the CD16 and NKG2D receptors are two of the most studied. In fact, a significant majority of NK cells in vivo are CD16-positive. Specifically, CD16 is an Fc c receptor found on the surface of NK cells, 7-10 which binds to the constant region of immunoglobulin G (IgG) antibodies. Due to its affinity for antibodies, CD16 is necessarily required for antibody-dependent cell-mediated cytotoxicity (ADCC), a mechanism for lysing target cells through antibodies. This feature of the CD16 receptor has been integral for designing bi-and tri-specific killer engagers (BiKEs and TriKEs), 11,12 which are engineered antibodies that traffic NK cells to target cells for cell killing. NKG2D belongs to the CD94/NKG2 family of receptors and has been found on NK cells as well as T cells. [13][14][15] Unlike CD16's ubiquity in ADCC, NKG2D is specific as it recognizes and binds to induced self-antigens [e.g., MHC class I polypeptide-related sequence A (MICA)] on the surface of cells. These antigens communicate to NK cells that the diseased cell should be lysed. This implicates NKG2D in the elimination of diseased cells, including cancer cells. Excitingly, NKG2D serves as a focal point for many lines of research in targeted therapies [15][16][17][18] due to its affinity for tumor-associated antigens.
While CD16 and NKG2D are activated under different biological scenarios, they activate a similar set of downstream signaling molecules [7][8][9][10]19 that mediate the secretion of GZMB and PRF1. Upon binding to their cognate ligands, antibodies, or antigens, CD16 and NKG2D promote activation of the Src family kinases (SFKs) through the intracellular adaptor molecules CD3f and DAP10, respectively. The activation of SFKs leads to the phosphorylation of downstream signaling species PLCc, Vav, SLP76, Akt, and Erk, as well as phosphatases SHP and SHIP. [8][9][10]20,21 The phosphatases inhibit the activation of the signaling intermediates and, thus, prevent cell activation. The activation of Vav and SLP76 is critical for actin remodeling and formation of the immunological synapse 4,8,22 between the NK cell and the target cell. Moreover, phosphorylation of Akt and Erk has been correlated with cell survival and proliferation, respectively. 7,23-25 Studies have shown 3,4,8,9,22,[26][27][28] a strong correlation between phosphorylation of the signaling molecules Vav and PLCc, secretion of GZMB and PRF1, and NK cell cytotoxicity, suggesting that the activation of these molecules precedes target cell death.
Several lines of research [29][30][31] have reported NK cell exhaustion as a consequence of over-stimulation of NK cell receptors. NK cell exhaustion is a phenotype characterized by a decrease in NK cell effector functions (e.g., GZMB/PRF1 secretion) 32 even with more receptor stimulation. In addition, NK cell exhaustion is correlated with a decrease in the density of stimulatory receptors. 29 Sanchez-Correa et al. 33 found a downregulation of DNAM-1 in primary NK cells exposed to DNAM-1 ligands CD112 and CD155 expressed by leukemic cells in vitro, leading to a dampened immune response upon subsequent stimulation of DNAM-1. Paul and colleagues 34 observed a reduced quantity of NKG2D-positive NK cells in the tumor microenvironment (TME) of B16F10-induced melanoma in C57BL/6 mice in contrast to NK cells in the periphery of the tumor site. Moreover, Paul et al. 34 discovered that NK cells within the TME had reduced levels of PRF1 mRNA, as well as lower levels of the cytokine IFN c and CD107 (a marker of degranulation), compared to NK cells in the periphery, suggesting that the TME can co-opt NK cells and promote the exhausted phenotype. In the clinical setting, a decrease in the quantity of NK cell stimulatory receptors and effector molecules has been correlated with poor prognoses for pancreatic, gastric, and colorectal cancer patients. 35,36 It follows that the efficacy of NK cell-based adoptive cell therapies will be limited unless NK cells can be modified to overcome the effects of exhaustion. However, there appear to be many mechanisms that induce NK cell exhaustion, 32 implying that a single approach may not work to prevent exhaustion. In addition, given the extensive cascade of signaling reactions that mediate NK cell degranulation and the complex interactions that influence NK cell exhaustion, it is not clear which strategies can be combined to reduce exhaustion or prevent it altogether.
Instead of preventing NK cell exhaustion, we could potentially implement strategies that promote its opposite effect, that is, the continuous secretion of cytolytic molecules. This approach may counterbalance the effects of exhaustion, causing the cell to be more robust to stimulation and, thereby, allowing the NK cell to kill more target cells. Still, it is not clear how to optimally achieve this objective due to the nonlinearities in cell signaling and activation. Excitingly, mathematical models are useful for providing quantitative insight into complex biological processes, including intracellular signaling leading to immune cell activation. 37 For example, mathematical models of NK cell signaling have contributed to our understanding of NK cell activation: the identification of Vav as the signaling molecule where stimulatory and inhibitory signals integrate to determine activation; 28 the explanation of how weak-affinity stimulatory receptors can counterintuitively inhibit NK cell activation in a non-monotonic manner; 38 the revelation that NK cell signaling occurs at a faster timescale than tumor cellkilling and how modifying antibody concentrations can bridge the two processes closer in time; 39 and the elucidation of strategies that increase the likelihood of NK cell activation 19 by amplifying the amount of phosphorylated signaling intermediates. Indeed, mathematical models can be applied to address a diversity of research questions, especially in determining optimal strategies, which can save experimental researchers' time and resources.
In this work, we have applied a mathematical model to investigate strategies that counteract NK cell exhaustion. We first modified our previous model of NK cell signaling 19 to include GZMB and PRF1 secretion. The model was first calibrated and validated using published data. 26 We then performed a global sensitivity analysis, revealing that activation of a particular phosphatase strongly influences NK cell secretion of cytolytic molecules. With that information in hand, we simulated the effects of reducing the phosphatase's impact, alone and in combination with increasing production of cytolytic molecules, for different rounds of stimulation. We found that the optimal strategy for maximizing secretion of cytolytic molecules depends on how many times the NK cells are stimulated: for fewer rounds of stimulation, inhibiting phosphatase activity leads to more secreted GZMB and PRF1; in contrast, for many rounds of stimulation, the production of the cytolytic molecules becomes essential. In conclusion, we constructed a mathematical framework describing the effector function of NK cells, which can aid researchers interested in engineering robust NK cells for clinical applications.

NK cell degranulation model can reproduce experimental observations
We generated a mathematical model of NK cell degranulation by incorporating the dynamics of GZMB and PRF1 into our previous model of NK cell signaling. 19 In total, the model consists of downstream signaling reactions from the receptors CD16 and NKG2D that ultimately mediate the secretion of the cytolytic molecules GZMB and PRF1. Specifically, the model consists of a system of nonlinear ordinary differential equations (ODEs) that predict the concentration of the receptors, signaling intermediates, and the cytolytic molecules; explicit equations, initial conditions, and parameters of the model are provided in supplementary material file S2. When the CD16 and NKG2D receptors are stimulated, they activate the cell via a cascade of reactions ( Fig. 1): activation of the Src family kinases (pSFK), facilitated by the ligand-bound phosphorylated receptors (pCD16 or pNKG2D), mediates the activation of the Akt, SLP76-Vav-Erk, and PLCc pathways. In particular, we assume that pVav and pPLCc mediate the secretion of GZMB and PRF1, as they have been correlated with NK cell cytotoxicity 29 and a discharge of intracellular calcium ions, 1,9,10,40,41 which mediate exocytosis, respectively. Although Akt and Erk are necessary for cell survival and proliferation, [8][9][10]41 respectively, their contribution to NK cell secretion of GZMB and PRF1 remains elusive. As it stands, given the available data, we consider pVav and pPLCc as the mediators for NK cell secretion.
The model was calibrated to data from the study by Srpan et al. 26 using a Bayesian perspective to parameter estimation; 47,48 namely, we implemented the Metropolis-Hastings (MH) algorithm (see the methods) to sample from the posterior distribution of the parameters conditional on the data. In brief, seven parameters (see Table I) were estimated 200 times using randomized initial guesses. We also tested the model predictions using a separate validation dataset. The combined error for each run can be found in Fig. S1. Since the marginal posterior distribution of each parameter for the best twenty runs was almost identical, we chose to simulate the model using the results from the best run (Run 1 in Fig. S1). The trace plots for each parameter in the best run are shown in Fig. S2, where the value of each parameter is plotted as a function of the iteration of the MH algorithm. This diagnostic of the parameter estimation shows that each parameter converges to a stationary distribution, albeit at different iterations. We simulated the model 1000 times using the final 1000 iterations of each parameter from the best run (Fig. S3) and compared the results with both the training and validation data. Excitingly, the model simulations are in good agreement with both the training and validation data (Fig. 2). The model can reproduce the time evolution of secreted PRF1 mediated by stimulation of NKG2D [ Fig. 2(a)] and CD16 [ Fig. 2(b)]. In addition, the model can replicate sequential stimulation data where NK cells are stimulated via one pathway for two consecutive rounds, each for 60 Fig. 2(k)], and these model simulations agree with experimental measurements as well. In addition, we demonstrate in Fig. S4 the model trajectories with respect to the experimental data and conditions.
To investigate the effects of using a different combination of training and validation data, we randomized the training and validation datasets by systematically switching the datasets from training to validation, and vice versa, to create sixteen different combinations including the original combination that we started with. For all combinations, we estimated the parameter values using the MH algorithm eight times using eight random initial guesses. Our results (Fig. S5) indicate that the best fit (Run 1 from Fig. S1) led to a smaller total error when compared to all combinations (including the original, which we re-estimated). Interestingly, combinations 4, 8, and 12 in Fig. S5 generated prediction errors similar in magnitude to the original combination. We found that the parameter distributions (from the final 1000 iterations of the MH algorithm) from combinations 4, 8, and 12 were comparable to the original combination with a significant overlap (Fig. S6). This suggests that the combinations that yield better fits to the data converge to similar parameter distributions during estimation, implying that the parameter estimates are robust with respect to variations in the model calibration process. Hence, we move forward with our original best fit and simulate the model using the parameter distributions found in Fig. S3, which overlap with the distributions in Fig. S6. Overall, we generated an experimentally validated mechanistic model of NK cell degranulation, which can recreate results from the study by Srpan et al. 26 in a range of different experimental conditions. We next apply the model to determine which features, when perturbed, robustly maximize the amount of secreted PRF1 and GZMB.
FIG. 1. Natural Killer cell signaling model schematic. Signal propagation flows from left to right, starting with the interaction between the ligand (Rtx or MICA) and the receptor (CD16 or NKG2D) that leads to the activation (phosphorylation) of the receptor-ligand complex (pCD16 or pNKG2D). Then, this complex can promote phosphorylation of the Src family kinases (pSFK), which activates the signaling intermediates as well as the inhibitory phosphatases (SHP and SHIP). Finally, the synergistic activation of pVav and pPLCc influences the release of granzyme B (GZMB) and perforin-1 (PRF1). Arrows indicate stimulation, whereas red crossbars signify inhibition.

Parameter
Units Description k deg min À1 Degradation of phospho-species k int CD16 min À1 Internalization and degradation rate of phospho-complex CD16 k int NKG2D min À1 Internalization and degradation rate of phospho-complex NKG2D k degran CD16 lM À2 Â min À1 Degranulation rate under CD16 activation k degran NKG2D lM À2 Â min À1 Degranulation rate under NKG2D activation k 0 min À1 Crosstalk from phospho-complex CD16 to NKG2D k 1 min À1 Crosstalk from phospho-complex NKG2D to CD16

Inhibition of pSHP maximizes PRF1 and GZMB secretion in silico
In order to understand which model parameters (inputs) leads to an increase in the predicted secretion of GZMB and PRF1 (outputs), we performed an entropy-based sensitivity analysis on the model (see the methods). This technique, as shown by L€ udtke et al., 50 measures the degree of mutual information shared between the model inputs and outputs. The greater the amount of shared mutual information between an input and output, the more sensitive the output is to that input. This metric is defined as the sensitivity index of the parameter (see the Methods). Briefly, we varied the model parameters 50% above and below their mean value and then drew 250 uniformly distributed samples to simulate the model and generate distributions for the amount of GZMB and PRF1 after 60 minutes of receptor stimulation. With these probability distributions for the model inputs and outputs, we are able to approximate the conditional entropies needed to compute the sensitivity indices.
Interestingly, only a select few parameters were shown to have a large sensitivity index in regard to GZMB and PRF1 secretion ( ] strongly influences the amount of GZMB and PRF1 secretion. In fact, the total effect of the most influential parameter, describing the catalytic rate of pSHP activation, shares over 70% of the information observed in GZMB and PRF1 secretion. Given the causal structure of the model, we can infer that blocking the pSFK ! pSHP reaction can yield more GZMB and PRF1 secretion since that would prevent the phosphatase from inhibiting signal transduction. Similar causal inferences can be made for other parameters in the subgraph in Fig. 3(c). Thus, we intervened on these parameters (individually) by modifying their values and simulating the model to measure the percent change in GZMB and PRF1 secretion from baseline compared to the modified case.
Since there are four variables of interest in this influential subnetwork (pSFK, pSHP, pVav, and pPLCc), each with a rate of activation and deactivation, we intervened on eight different parameters to enhance or inhibit the phosphorylated species. Specifically, the catalytic rate constants were varied from baseline. Then, the model was subsequently simulated to quantify their impact on the percent change in GZMB and PRF1 secretion after 60 min of receptor stimulation  . These results suggest that disrupting the incoherent feedforward loop (IFFL) pSFK ! pSHP a pX, where X is either Vav or PLCc, leads to more GZMB and PRF1 secretion. In summary, we interrogated the model and found that the amount of cytolytic molecule secretion is strongly affected by the pSFK ! pSHP edge in the subgraph depicted in Fig. 3(c). The optimal synNotch system depends on the number of rounds of NK cell stimulation The above results indicate that a strategy to increase GZMB or PRF1 secretion is to inhibit pSHP activation. While this does lead to more GZMB and PRF1 secretion in silico, it almost completely depletes the intracellular pool of cytolytic molecules (Fig. S9). Depleting intracellular pools of GZMB and PRF1 makes the NK cell less likely to secrete these cytolytic molecules upon subsequent stimulation since the timescale over which the pool is replenished (i.e., protein production) is longer than the timescale over which the cell would be stimulated. 39,65 Ultimately, this causes NK cells to become less cytotoxic over time. Therefore, to maximize GZMB and PRF1 secretion over multiple rounds of stimulation, we need to induce the production of these molecules, in addition to reducing SHP activation. Fortunately, the field of synthetic biology 52-54 provides tools that can be applied to address this issue. Multi-cistronic plasmids, which express two or more genes, can be used to promote expression of both the GZMB and PRF1 genes in NK cells and, thereby, replenish the intracellular pool of the cytolytic molecules. When coupled with the inhibition of pSHP, this strategy may increase both cytolytic molecule secretion and production, enabling NK cells to continuously secrete cytolytic molecules over multiple rounds of stimulation.
Excitingly, the synthetic Notch (synNotch) signaling system 55 can be applied to simultaneously (1) inhibit pSHP and (2) increase the intracellular pools of the cytolytic molecules when the NK cell senses a danger signal (e.g., MICA). In brief, the synNotch system can be constructed via genetic modifications. Once the cell expresses the synNotch receptor, it is able to promote a specific function 52-54,77 that can be tailored to a specific stimulus [ Fig. 4(a)]. The NK cell signaling network promoted by the endogenous receptors CD16 and NKG2D combined with the synthetic Notch system is illustrated in Fig. 4(b). We apply this augmented modeling framework as a mechanism for controlling the production and secretion of the cytolytic molecules (i.e., the specific functions) when the NK cell binds to CD16 and NKG2D ligands (i.e., the specific stimuli). For inhibiting pSHP, we simulate the expression of long non-coding RNA (lncRNA) targeting SHP, which impedes SHP from binding to its targets, thereby reducing its ability to inhibit signaling 58,59,61 (see Methods). At the same time, we consider increasing the expression of GZMB and PRF1 using the synNotch system. Collectively, this approach can grant NK cells the ability to maximally secrete cytolytic molecules upon repeated receptor stimulation.
We study the case where the extracellular domain of the synNotch receptor has the same binding kinetics as the endogenous receptor (see Table II). In this case, the synthetic and endogenous pathways are independent yet complementary to one another: the synthetic path produces the cytolytic molecules, while the endogenous pathway secretes them. Although this strategy seems intuitive, the optimal amount of synNotch receptor is difficult to deduce a priori given the competing dynamics for the ligand. Moreover, it is unclear how much of each plasmid is required to maximally secrete GZMB and PRF1 when considering the number of rounds of stimulation. Since transcription and translation are demanding biological processes, 51,54,76 which are required to express the synNotch receptor in addition to the lncRNA and cytolytic molecules, we consider using the absolute minimal amount of exogenous material needed to achieve maximum secretion. Taken together, we aim to find the optimal levels of lncRNA-coding plasmid targeting SHP (lncSHP), cytolytic molecule-coding plasmid, and the synNotch receptor needed to maximize GZMB and PRF1 secretion over many rounds of stimulation while using the absolute minimal amount of synthetic material (see the Methods).
Interestingly, the optimized synNotch system has different characteristics when paired with CD16 [ Fig , it is optimal to do nothing. That is, the amount of effort required to express the synNotch system outweighs the gain in performance (i.e., increased secretion). This is due to the fact that the time delay inherent in transcription and translation. That is, two rounds ($2 h) of receptor stimulation is not enough time to appreciably change the concentration of intracellular GZMB and PRF1 via the synNotch system. This holds true for NKG2D as well [ Fig. 4(f)]. In contrast, when CD16 is stimulated for 3-5 rounds, the optimal amounts for each plasmid and R 0 (initial value of synNotch) are at their maximum values. At this stage, the synNotch system is fully operational. Surprisingly, as we continue to increase the number of rounds of CD16 stimulation, we found that the optimal amount of lncSHP-coding plasmid decreases precipitously by several orders of magnitude, while the optimal amounts of the cytolytic molecule-coding plasmid and R 0 remain at their maximum values. Intriguingly, the model prediction reveals that inhibition of SHP is not optimal in the long-term, as the optimal amount of lncSHP plasmid is effectively zero. This is because too much inhibition of SHP leads to an accumulation of phospho-proteins, which increases The ligands (red) can bind to the synNotch receptor, which allows the complex to change in conformation to allow constitutively expressed, membrane-bound proteases (green rectangle) to cleave the peptide link between the synNotch receptor and the transcription factor (TF, blue). This allows the TF to bind to its binding site (cyan) on the plasmid, which initiates gene (orange) transcription. (b) Schematic of interaction between the endogenous and synthetic pathways. The arrows and crossbars represent stimulation and inhibition, respectively. The objective function is minimized over various rounds of stimulation of (c) CD16 and (f) NKG2D. For each round of stimulation, the optimal amounts of lncSHP plasmid (red markers), cytolytic molecule plasmid (blue markers), and synNotch receptor (yellow markers) are shown. (d) and (e) The performance from CD16 stimulation with and without controllers. The secretion of (d) GZMB and (e) PRF1 via the CD16-synNotch pair pathway. (g) and (h) The performance from NKG2D stimulation with and without controllers. The secretion of (g) GZMB and (h) PRF1 via the NKG2D-synNotch pair pathway. No control, optimal values from round 1 for CD16 and NKG2D (gray squares); max control, optimal values from round 3 for CD16 and round 11 for NKG2D (green squares); semi-max control, optimal values from round 6 for CD16 and round 15 for NKG2D (orange squares). Markers, mean model prediction; error bars, one standard deviation.
the velocity of phospho-protein decay 10,78,79 (see Table I, k deg ). Moreover, since there is no synthesis reaction of the inactive signaling species in our model, the degraded phospho-species reduces the total amount of available molecules for the next round of signaling. Thus, for many rounds of stimulation, the phosphatase counterintuitively is needed to maintain the availability of proteins for subsequent signaling.
The optimal synNotch system when coupled with NKG2D has subtle yet important differences. The model predicts that when NKG2D is stimulated for three rounds or more, the optimal amount of cytolytic molecule-coding plasmid is at its maximum [ Fig. 4 Fig. 4(c)]. Unlike CD16, however, the optimal amount of lncSHP-coding plasmid remains at a maximum value from 3-14 rounds of stimulation, and the optimal concentration of R 0 does not immediately reach its maximum value. This difference between the optimal amounts of R 0 is due to the initial amounts of CD16 and NKG2D receptors (38 and 0.3 lM, respectively), compared to the maximum value that R 0 can have, 10 lM 76 (see the Methods section). We interrogated the model and found that the concentrations of CD16 and NKG2D lead to this distinction. The high concentration of the endogenous CD16 receptor compared to the concentration of the synNotch receptor means that CD16 will outcompete the synNotch receptor. This helps clarify why the optimal value of R 0 is the maximal value it can take on. In comparison, for NKG2D stimulation, where the initial amount of endogenous receptor is less than the maximal amount of R 0 , a large R 0 would lead to less secretion by directing the input signal more toward cytolytic molecule production. Moreover, since the optimal R 0 is small when considering 10 or fewer rounds of NKG2D stimulation, the indirect impact of SHP inhibition on phospho-protein decay is also small, thus allowing the amount of lncSHP-coding plasmid to remain at the maximum for more rounds of NKG2D stimulation. Nevertheless, as the rounds of NKG2D stimulation increase, where the optimal R 0 approaches its maximal value, the optimal amount of lncSHP-coding plasmid is reduced by several orders of magnitude -similar to the case for CD16 stimulation. In summary, we found that the optimal synNotch system not only depends on the number of rounds of receptor stimulation but also is different when considering whether the CD16 or NKG2D pathway is being stimulated.
The predicted optimal CD16-and NKG2D-synNotch pairs show qualitative differences in cytolytic molecule secretion Given the optimized synNotch system, we next simulated the model to predict how the secreted amount of GZMB and PRF1 changes, compared to the baseline case without the synNotch receptor. We simulated the model using the optimal amounts of each plasmid and the initial synNotch receptor from different rounds of stimulation to observe any differences in GZMB and PRF1 secretion. The results of our optimization analysis show three separate clusters of optimal conditions, which depend on the number of rounds of receptor stimulation: (1) no synNotch system, (2) maximal amount of both plasmids and synNotch, and (3) no lncSHP-coding plasmid but maximal amount of cytolytic molecule-coding plasmid and synNotch. For simplicity, we label these clusters as "no control," "max control," and "semi-max control," respectively. For CD16, we used the optimal amounts from Fig. 4(c) for rounds 1, 3, and 6, respectively, for these three conditions; whereas for NKG2D, we used the optimal amounts from Fig. 4(f) for rounds 1, 11, and 15, respectively.
When considering GZMB secretion [ Fig. 4(d)], the model analysis shows that using the amounts of synNotch receptor and plasmids optimized for one round of stimulation (in this case, none at all), more GZMB is released compared to using the amounts of receptor and plasmids optimized for multiple rounds of stimulation in the first few rounds of stimulation. Specifically, not having the synNotch receptor is best for up to one round of stimulation of the CD16 pathway and up to two rounds of stimulation of the NKG2D pathway [ Fig. 4(g)]. This counterintuitive result that no control is optimal is due to the competition for ligand between the endogenous and the synNotch  Interestingly, when the synNotch system is absent, we found that one round of NKG2D stimulation leads to much more secreted PRF1 than CD16 stimulation (3.8 vs 2.0 pmol, respectively), corroborating the findings in the study by Srpan et al.. 26 However, for the second round of stimulation, the amount of secreted PRF1 reduces significantly for NKG2D stimulation to 1.4 pmol, while CD16 stimulation yields 1.5 pmol of secreted PRF1. For the subsequent round of stimulation, NKG2D yields 0.8 pmol, whereas CD16 produces 1.5 pmol. This initial large burst in NKG2D-mediated PRF1 secretion can be explained by the estimated degranulation parameter k degran NKG2D (Fig. S3), which is approximately three times larger than k degran CD16 , and thus, the velocity of secretion is faster via the NKG2D path. Taken together, in the absence of a synthetic pathway, NKG2D leads to a large but transient secretion of PRF1, whereas CD16 produces a small but steady secretion of PRF1 after two rounds of stimulation in silico.
In contrast, when we consider more rounds of stimulation, it is best to use the synNotch system. Using the max control set The model again predicts differences between the CD16-synNotch pair and the NKG2D-synNotch pair. When using the max control set, the secreted amount of GZMB mediated by the CD16-synNotch pair decreases to zero after 7 rounds of stimulation [Fig. 4(d)]. The NKG2D-synNotch pair, on the other hand, continues to secrete GZMB up to 14 rounds of stimulation. The model predicts similar differences when considering secreted PRF1, where secretion of PRF1 promoted by stimulation of CD16 goes to zero after 8 rounds of stimulation. In contrast, PRF1 secretion promoted by NKG2D stimulation goes to zero after 16 rounds of stimulation. The differences in the concentration between synNotch and the endogenous receptors are mainly responsible for the differences in the responses. Given that R 0 is greater than the initial value of NKG2D (0.3 lM), more of the signal is veered toward cytolytic molecule production due to competition for the ligand. This explains NKG2D-synNotch's sustained response to stimulation. The initial concentration of CD16 (38 lM), in contrast, is greater than R 0 , and therefore, more of the signal is shifted to cytolytic molecule secretion, resulting in a large but transient response. Taken together, in the presence of synNotch, the qualitative features of CD16-vs NKG2D-mediated secretion of the cytolytic molecules changed compared to the no control case: the CD16 path now produces a large but transient response to stimulation, while the NKG2D pathway yields a small but sustained response to stimulation.
Interestingly, the optimal synNotch system does not indefinitely increase the amount of secreted cytolytic molecules as we increase the frequency of NK cell stimulation. In fact, when the number of rounds of receptor stimulation becomes large [Figs. 4(d), 4(e), 4(g), and 4(h)], the secreted amount of cytolytic molecules decreases. We interrogated the model to better understand this prediction. While the intracellular pool of GZMB [ Fig. S10(a)] and PRF1 [ Fig. S10(b)] continues to increase as we increase the number of rounds of stimulation, the intracellular amount of pPLCc [ Fig. S10(c)] and pVav [ Fig. S10(d)] eventually decreases to zero. Since the secretion of GZMB and PRF1 is mediated by pPLCc and pVav in our model, it follows that once the amount of pPLCc and pVav becomes negligible, so too does the secretion of GZMB and PRF1. Thus, at higher frequencies of stimulation, the depletion of intracellular signaling species is predicted to be a limiting factor in cytolytic molecule secretion.
To assess the robustness of the above results, we performed a sensitivity analysis on the optimized model that includes the synNotch pathway (see Table II). Similar to our previous sensitivity analysis, we varied each parameter by 50% above and below its baseline value to create a uniform distribution. Next, we drew 250 samples and simulated the model to determine how sensitive the secretion of GZMB and PRF1 is to the parameters under stimulation of CD16 (Fig. S11) or NKG2D (Fig. S12). This sensitivity analysis was performed using the max control set (3 and 11 rounds of stimulation for CD16 and NKG2D, respectively) and the semi-max control set (6 and 15 rounds of stimulation for CD16 and NKG2D, respectively) of optimal results from Figs. 4(c) and 4(f). For CD16, we found that the total order sensitivity indices of each of the parameters in the synNotch system, for GZMB secretion from the optimal results for the max control set (round 3), can at most explain 32% of the information in the output [ Fig. S11(a)]. Moreover, when compared to the parameters in the endogenous pathway [ Fig. S11(a); orange vs black], the synNotch parameters are much less relevant, meaning that the secretion of GZMB is robust to our assumptions related to the synNotch signaling model. The parameters that are found to be the most influential are the same as those found in Fig. 3(a), which belong to the endogenous pathway. Similarly, these conclusions hold for the semi-max control case of GZMB secretion [ Fig. S11(c)] as well as for PRF1 secretion [Figs. S11(b) an S11(d)] and GZMB and PRF1 secretion promoted by the NKG2D pathway (Fig. S12). To further verify these results, we intervened on the R 0 and k on parameters in the synNotch pathway (which are the most influential parameters from the synNotch system) as we presented above [ Fig. 3(d)] and compared the results with the decreasing activation rate of pSHP (Fig. S13). We found no significant change in the secretion of GZMB or PRF1. In conclusion, we augmented our baseline model with a synNotch signaling system and found that the secretion of the cytolytic molecules can be enhanced with this added system, and the predicted impact of the synNotch system is not significantly affected by the parameters and assumptions we implemented.
Transcription factor cooperativity impacts the performance and optimality of the synNotch system Since transcription and translation restrict the rate of cytolytic molecule production, we investigate in silico the optimality of the synNotch system as we modify the transcriptional dynamics. The results from Figs. 4(c) and 4(f) show that the synthetic circuit is not needed if NK cells are stimulated for less than three hours due to the time delay in transcription and translation. Studies have demonstrated that modulating transcription 80,81 can enhance the rate of protein production. With this knowledge at hand, we alter key parameters that may increase RNA production. In this work, we represent the rate of transcription of DNA as a Hill function of the transcription factor (TF) (see Table II where u is the amount of plasmid and k trpt is the transcription rate. By taking u and k trpt to be fixed, we observe that the rate of RNA production increases as we decrease the parameters n (the Hill coefficient) and K M (the affinity between the DNA molecules and the TF). We consider the case of varying the Hill coefficient only, although a similar argument can be made for the affinity parameter. Briefly, the Hill coefficient captures the sharpness of a response (i.e., RNA production) to a given input [82][83][84] (i.e., the TF). It is possible to change the rate of RNA production by varying n. For example, decreasing n should increase the velocity of RNA production. Indeed, if we take, for example, n 0 ¼ def 0:5n and n 1 ¼ def 1:5n, then, where r t ð Þ is the rate of RNA production with n ¼ 1, and similarly, for all time t. It follows that if we decrease n, we should expect to find an increase in performance of the synNotch system. To investigate the effects of the Hill coefficient, we repeat the analyses and simulations as we did above for two new cases: (1) n ¼ 0.5 and (2) n ¼ 1.5 to represent "negative" and "positive" cooperativity, respectively. As expected, we found that when n < 1, we enhance the performance of synNotch since the optimal amounts of lncSHP-coding plasmid [ Fig. 5(a)] and cytolytic molecule-coding plasmid [ Fig. 5(b)] are at their maximal values and R 0 [Fig. 5(c)] is near-maximal even with one round (i.e., 1 h) of CD16 stimulation. By increasing RNA production, we increase the performance term in the objective function relative to the effort term, and therefore, the synthetic circuit will be useful in short-term receptor stimulation. Contrastingly, when n > 1, we reduce the necessity of synNotch as evidenced by a delay in intervention [compare the lightest shade (n < 1) with darker shade (n > 1) in Figs. 5(a)-5(c)]. These results hold for the NKG2D pathway as well [Figs. S14(a)-S14(c)]. However, the data suggest that augmenting the velocity of RNA production only affects the optimality of synNotch when considering five rounds of stimulation or fewer in the case of CD16. Note that in Figs. 5(a)-5(c), when the number of rounds of CD16 stimulation increases, the optimal amounts of synNotch components are nearly identical for the different values of n. In summary, the optimality of synNotch is sensitive to transcriptional dynamics. Given the optimal synNotch components for the different levels of TF cooperativity, we simulate the model to predict and compare the amount of secreted cytolytic molecules in the different cases. We simulate the model using the control sets (i.e., no control, max control, and semi-max control) we defined previously. Here, for the CD16 pathway, the no control, max control, and semi-max controls corresponded to using the optimal synNotch components from rounds 1, 3, and 6 from Fig. 4(c). We apply the model by using the predicted optimal values from the same rounds with the different Hill coefficients. When we consider GZMB secretion [Figs. 5(d)-5(f)], we found that when n ¼ 0.5 (lightest shade), the performance of synNotch is superior to all other cases, irrespective of the control set. These results hold true for PRF1 secretion [Figs. 5(g)-5(i)] and the NKG2D pathway [ Fig. S14(d)-S14(i)]. Interestingly, in the case of NKG2D, we found that the qualitative nature of the synNotch-NKG2D pair resembles that of the synNotch-CD16 pair when n ¼ 0.5 [compare the lightest shade in Figs. S14(e) and S14(f) with Figs. 5(e) and 5(f)], suggesting that increasing RNA production shifts the response from small-butsustained to large-and-transient [compare the lightest shade with intermediate shade in Figs. S14(e), S14(f), S14(h), and S14(i)]. In conclusion, we found that TF negative cooperativity is predicted to augment the performance of synNotch, whereas positive cooperativity is predicted to reduce it.

The synthetic endogenous pathway pairs promote a biphasic steady state response to stimulus
The analysis of dose-response relationships is helpful in understanding the link between the input and output of a system. Here, we apply the model, with and without the controllers, to determine how the steady state amount of secreted GZMB and PRF1 (i.e., response) varies as we vary the amount of ligand (i.e., dose). More specifically, we vary the concentration of Rituximab (Rtx) and MICA from 10 À5 to 10 lM and simulate the model for one round of stimulation until the steady state is attained. That is, instead of stimulating the receptors for k rounds (60 min each) with kÀ 1 intermediary washing steps (15 min each), we simulate the stimulation of the receptors for one round but for a sufficiently long time-horizon (i.e., > 40 h). Given the signaling dynamics, the only species with non-trivial steady state values are the cytolytic molecules; the amount of each intracellular phospho-species goes to zero as time approaches infinity as there is only a sink in the current model represented by the degradation term (see k deg , Table I). We simulate the steady state response with the base case where the Hill coefficient (n) is equal to one.
As expected, the steady state response of secreted GZMB via the CD16-synNotch pair [ Fig. 6(a)] is larger when the controllers are present (compare green and orange curves with the gray curve). Interestingly though, the response is biphasic with respect to the magnitude of stimulus; that is, as we continue to increase the amount of Rtx beyond 0.7 or 0.2 lM-in the max and semi-max control settings, respectively-the steady state levels of GZMB do not increase but, in fact, decrease. Given the endogenous and synthetic signaling dynamics, the model predicts the existence of a range where synNotch is useful. In the case of GZMB [ Fig. 6(a)], the benefits of synNotch are observed when Rtx is between 10 À3 and 10 lM and peak at approximately 0.7 and 0.2 lM for the max and semi-max control cases, respectively. The data suggest that a stimulus level too small or too large nullifies the effects from synNotch. We observe a similar characteristic response when considering the steady state levels of PRF1 via the CD16-synNotch pair [ Fig. 6(b)] or via the NKG2D-synNotch pair [Figs. 6(c) and 6(d)].
Next, we probed the model to better understand the cause of the biphasic, steady state response. In Figs. S15(a)-S15(d), we demonstrate the time evolution of secreted GZMB via the CD16-synNotch pair at four distinct levels of Rtx: 10 À2 , 10 À1 , 1, and 10 lM, respectively. Also, we show the time evolution of secreted PRF1 in Figs. S15(e)-S15(h). We use the same control sets as we defined above. In these figures, we observe that when Rtx is 10 À1 lM, the steady state response is largest for both species when the controllers are present. Intriguingly, as we increase the magnitude of Rtx, note that the time required to reach the steady state decreases [compare, for example, Figs. S15(a) with S15(b)-S15(d)]. More importantly, phospho-CD16-the species that mediates the activation of the endogenous pathway, thereby leading to cytolytic molecule secretion-continues to signal well over 20 h when Rtx is equal to 10 À2 lM [ Fig. S15(i)] and close to 8 h when Rtx is equal to 10 À1 lM [ Fig. S15(J)]. However, once the concentration of Rtx is 1 or 10 lM (Figs. S15(k)-S15(l)], the concentration of phospho-CD16 undergoes a large rapid increase but decreases precipitously after a 3or 1-hour stimulation period, respectively. This response is attributable to the signaling dynamics regulating phopsho-CD16, where the velocity of internalization and ligand recycling (see k int CD16 , Table I) is proportional to the amount of phospho-CD16. That is, degradation of phospho-CD16 / k int CD16 pCD16 t ð Þ, where k int CD16 is between 1.4 and 1.5 min À1 [ Fig. S3(b)]. Although a large concentration of Rtx shortens the time to reach the steady state, it reduces the magnitude of secreted GZMB and PRF1 at steady state because the amount of phospho-CD16 depletes faster and, therefore, ceases to promote the signal for secretion of cytolytic molecules. These conclusions are shown to be consistent with the NKG2D-synNotch pair as well (Fig.  S16). In summary, our results indicate that the long-term benefits of synNotch are biphasic with respect to the amount of stimulus due to depletion of the phosphorylated endogenous receptor.

DISCUSSION
In the present work, we constructed a mathematical model of NK cell secretion of GZMB and PRF1, which replicated experimental observations from the study by Srpan et al.. 26 The model was simulated to better understand which subnetworks strongly regulate the secretion of the cytolytic molecules as well as strategies for maximizing their secretion in silico. Furthermore, we investigated in silico how to enhance secretion of the cytolytic molecules over time. Specifically, we simulated the effects of engineering the NK cell to express the synNotch receptor that simultaneously enables production of GZMB and PRF1 and expression of an inhibitor of the phosphatase that significantly affects GZMB and PRF1 secretion. Our modeling revealed the following: as the number of rounds of stimulation increases, the optimal initial value of the synNotch receptor increases proportionally. This implies that if the NK cell is to be stimulated for multiple rounds, it is optimal to promote the synthetic pathway in order to increase production of GZMB and PRF1. In general, the production of the cytolytic molecules should be induced as maximally as possible. The optimal amount of lncSHP-coding plasmid, however, should switch from its maximum value to the minimum (0 copies per cell) as the number of rounds of stimulation increases, suggesting that SHP inhibition is only optimal in the short-term.
Moreover, negative cooperativity in transcriptional dynamics enhances RNA production and, thereby, augments the performance of synNotch. The long-term behavior (steady-state GZMB or PRF1 secreted) in response to varying ligand concentrations of the endogenous signaling pathway paired with the synthetic system is biphasic, implying that too small or too large of an input leads to a suboptimal performance. In total, our work presents a theoretical framework that can be used by researchers for engineering NK cells with applications to cell-based therapies.
The baseline model predicts that cytolytic molecule secretion is strongly dependent on the IFFL pSFK ! pSHP a pX, where X is either Vav or PLCc. Interestingly, this motif is not uncommon in biological networks [85][86][87][88] given its significant role in controlling cell activation by regulating signal transduction. Indeed, occluding the above subnetwork disinhibits signal transduction, allowing the NK cell receptor to continue signaling to its downstream mediators and, thereby, generate a greater response. In fact, there are many examples of increased cell activation using pharmacological inhibitors of phosphatases. 20,[89][90][91][92] Here, we observed a similar result by decreasing the catalytic rate constant that regulates the rate of SHP activation in silico. Still, given the complexity of the signaling network, the efficacy of inhibiting the phosphatase is not immediately obvious. Our global sensitivity analysis revealed the importance of the phosphatase, and simulating the effects of reducing phosphatase activity confirmed its impact, demonstrating the utility of the model. The synNotch signaling system is a powerful tool for engineering cells with novel capabilities. One advantage of this approach is that it grants the researcher the ability to program new pathways in cells using established methods in molecular biology; however, arriving at the optimal signaling circuit via experimentation alone can be cumbersome, expensive, and time consuming. In addition, it is often not clear if optimality was obtained. Excitingly, mathematical models of cell signaling systems can help address the above barriers. Namely, the question of optimality can be addressed when a model is applied to perform optimal control, given a set of controllers and an objective function. In the work presented here, the initial amounts of plasmidencoding lncRNA for the phosphatase, the plasmid-encoding cytolytic molecules GZMB and PRF1, and the synNotch receptor are analogous to controllers as they help steer the NK cell signaling network to secrete more cytolytic molecules. Once a solution is found to the optimization problem, we gain insight on how to optimally control NK cell degranulation.
Our findings demonstrate that the optimal strategy for maximizing the secretion of the cytolytic molecules is dependent on the number of rounds of stimulation that the cell will experience. We found that the optimal plasmid amounts have a bang-bang characteristic, meaning that they are either applied at maximal dose or not at all. This is a known characteristic for bound controls that appear linearly in the objective function. 69 In comparison, the optimal initial value of the synNotch receptor is more tunable and generally increases proportionally with the number of rounds of stimulation. Given that the synNotch receptor competes with the endogenous receptor for ligands, the increase in its initial amount reflects how much of the signal should be shifted to cytolytic molecule production and SHP inhibition vs secretion. Transcriptional dynamics are influential in determining the performance of the synthetic pathway, where decreasing TF cooperativity leads to an increase in the velocity of RNA production and, thus, an increase in cytolytic molecule production. Although engineering negative cooperativity may be difficult in practice, an analogous strategy would be to increase the affinity 93-95 between the TF and its DNA binding site on the plasmid (i.e., decreasing K M ).
The steady state performance of the synNotch system is sensitive to the amount of ligand. Specifically, we found that when the magnitude of stimulus is less than 10 À3 lM or greater than 10 lM, the added benefits of synNotch are nullified. In fact, the speed of degradation of the intracellular signaling molecules is proportional to the amount of stimulus, implying that a sufficiently large stimulus is suboptimal for cell signaling since it leads to rapid depletion of the molecules and, therefore, impedes cytolytic molecule secretion. On the other hand, when the amount of ligand is scarce-as in the case of tumor cells with low immunogenicity-endogenous and synthetic signaling will be negligible, and other strategies for enhancing secretion of cytolytic molecules should be investigated in these cases.
Our modeling and analysis predict the optimal amounts of controllers for specific rounds of stimulation. As it stands, however, it is not completely understood how many target cells an individual NK cell can kill in its lifespan; that is, how many times an NK cell would be stimulated. Prager et al. 3 recently showed that a given NK cell in vitro can kill up to six HeLa cells when confined to a microfluidic device. Elsewhere, others have demonstrated that NK cells are capable of serial killing, 3,26,[96][97][98] implying that they are able to undergo multiple rounds of stimulation. It remains to be seen if these results hold in vivo as well. While there is uncertainty in the exact killing potential, our analysis predicts that the highest secretion of GZMB and PRF1 occurs when we expect that NK cells will undergo many rounds of stimulation. The optimal synthetic biology solution is to have the cytolytic molecule-coding plasmid given at the maximal level, with high affinity between the TF and its DNA binding domain on the plasmid, in combination with a synNotch receptor to generate NK cells that continuously secrete effector molecules, even up to almost 10 rounds of stimulation. Although the precise numerical values of the optimal sets are subject to vary in practice, the qualitative predictions from the model are particularly useful, that is, the effects of using the maximum (or minimum) values of the plasmids and synNotch. Also, given that the intracellular pool of signaling species can limit the amount of secreted cytolytic molecules and that the production of each species may not be feasible in practice, the optimal strategies proposed here complement ongoing research aimed at increasing the population of NK cells or their proliferative capacity along with the strategies presented here. Taken together, these approaches may improve the efficacy of NK cell-based therapies.
We acknowledge that the model predictions are sensitive to the assumptions made on the model (see the Methods). We assume that an increase in the secretion of the cytolytic molecules will increase the likelihood of target cell death since these effector molecules are the mediators for apoptosis. Given the model was fit to data from the study by Srpan et al., 26 where the ligands were immobilized in each well of the 96-well plate, we assumed that the total amount of ligand in the system was fixed. This, however, need not be true in vivo where the ligands are subject to degradation, shedding, and clearance. The assumptions on the synNotch signaling pathway (see Table II) were implemented to simplify the optimal control analysis. Certainly, the extracellular domain of synNotch receptor does not need to have the same properties as the endogenous receptor; the binding and internalization kinetics can be different in practice. Furthermore, the affinity between the transcription factor and the plasmid can vary depending on what specific nucleotide sequences are used in the promoter site and the specific transcription factor. In the future, we can simulate the model where one plasmid has a higher affinity to the transcription factor than another plasmid, which may affect the optimal strategy. We also assumed that lncRNA can bind to both SHP and pSHP with equal affinity. Fortunately, these assumptions can be addressed by experimentation and parameter estimation to improve the precision of the model. Finally, we acknowledge that the solution to the optimization problem is sensitive not only to the model but also to the objective function. In particular, the constant q determined how much emphasis we placed on minimizing the given amount of exogenous material vs maximizing the cumulative amount of secreted cytolytic molecules. In our analysis, we placed an equal weight on both. Future research can address this issue by solving the optimization problem for a variety of desired outcomes (e.g., more emphasis on performance). We note that the data used for fitting were from a single published study. We acknowledge that using data from different sources would augment the training and validation of the model. However, there is a lack of studies involving repeated stimulation of NK cells where the number of ligands is quantified and standardized before cell stimulation, an input needed for model simulation. The model can be further validated as additional data become available.
Despite the limitations, the simulated results provide insight into NK cell secretion of cytolytic molecules GZMB and PRF1 mediated by the CD16 and NKG2D signaling pathways. We trained and validated a mathematical model by estimating the posterior distribution of the model parameters using the Metropolis-Hastings algorithm. An information-theoretic sensitivity analysis was subsequently performed on the model, from which we identified that the inhibition of SHP strongly influences the secretion of the cytolytic molecules. By incorporating a synNotch signaling system, we found that the optimal conditions for maximizing degranulation are dependent on the number of rounds of receptor stimulation: we found that SHP inhibition is not optimal in the long-term, while the production of the cytolytic molecules should be maximally induced almost all the time. Additionally, we found that negative cooperativity in transcription improves the performance in synNotch signaling and that the steady state response of the endogenous pathway combined with the synthetic system is biphasic with respect to the amount of stimulus. In conclusion, the current work provides actionable insight into engineering robust NK cells with applications to immunotherapies.

METHODS
We adopted the mathematical model of NK cell signaling from our previously published work, 19 which includes a system of nonlinear ODEs that describes the dynamics of receptors CD16, NKG2D, and 2B4, as well as their downstream signaling intermediates. The system of ODEs is integrated using the MATLAB function ode15s. We apply a range of computational methods (described below), and ethics approval is not required.

Construction of the NK cell degranulation model
Here, we focus on CD16 and NKG2D signaling, given extensive experimental data quantifying their roles in promoting cytolytic molecule degranulation. Furthermore, we expanded the model by incorporating the dynamics of GZMB and PRF1 to study and simulate strategies that optimize their secretion. This baseline model is provided in supplementary material file S1, and the list of model species, reactions, and parameters is provided in supplementary material file S2. A simplified representation of the baseline model is shown in Fig. 1, where connections between species represent reaction velocities and are governed by established Michaelis-Menten kinetics.
The baseline model contains 89 parameters and 40 species. Briefly, each receptor binds to its ligand and forms a receptor-ligand complex that can then become phosphorylated by basally active Src family kinases (SFKs). Then, the ligand-bound phosphorylated receptor (pC in Fig. 1) serves as the catalyst for converting SFK from a basally active state to a fully active state (pSFK). 29 Next, pSFK mediates the phosphorylation of PLCc, Vav, SLP76, Akt, and Erk, in addition to the phosphatases SHP and SHIP. 7,10,[40][41][42] The phosphatases SHP and SHIP provide negative feedback to the stimulatory network by dephosphorylating the phosphorylated signaling intermediates, including pC and pSFK. The initial concentrations of the model were taken from the literature, 43 and the kinetic parameters regulating the rate of phosphorylation and dephosphorylation reactions are taken from our published model. 19 In particular, the initial concentrations are assumed to be steady state values since the primary NK cells were pre-incubated in cell culture media, which excluded stimulatory ligands for CD16 or NKG2D prior to the measurement via mass spectrometry. 43 In addition to the model given by Makaryan and Finley, 19 we included a degranulation parameter for each pathway (CD16 and NKG2D) to account for actin remodeling, trafficking, docking, and exocytosis of the cytolytic molecules. 4,8,10,40,44 Given pVav's correlation with NK cell cytotoxicity 29 and pPLCc's role in the release of intracellular calcium ions, which are needed for exocytosis, 1,2,4,40-42 we used pVav and pPLCc as catalysts for GZMB and PRF1 secretion. Also, data from the study by Srpan et al. 26 demonstrated crosstalk between CD16 and NKG2D; specifically, stimulation of CD16 induced a slight increase in the amount of NKG2D, whereas stimulation of NKG2D slightly decreased the concentration of CD16. We included reactions to account for the crosstalk, based on either mass action kinetics or a nonlinear Hill equation. We also considered alternate models that included synthesis and decay reactions for the inactive signaling species since the data in the study by Srpan et al. 26 stimulated NK cells over a long timescale. Excitingly, all candidate models demonstrated a good agreement with experimental observations. We compared each model using the Akaike information criterion (AIC) that measures the quality of models given a dataset. We found a simpler model, which (1) had linear equations for the crosstalk reactions and (2) excluded synthesis and decay reactions for the inactive species, to be the preferred model (see Table III). Specifically, we calculated the relative likelihood of the models: exp 1 where i is the index number of the model. This value represents how likely or probable the model candidates are to the one with the lowest AIC score.

Data collection and processing
The mathematical model was trained using experimental data from the study by Srpan et al., 26 who used (1) primary NK cells in their experimental studies, (2) Rituximab as the ligand for stimulating CD16, and (3) MICA for the stimulation of NKG2D. The freely available online software WebPlotDigitizer (https://automeris.io/ WebPlotDigitizer) was used to extract mean and standard deviation data from plots shown in the study by Srpan et al. 26 In that study, the researchers stimulated primary NK cells in a 96-well plate using immobilized Rituximab and MICA for activation of CD16 and NKG2D, respectively. In addition, each well was coated with anti-PRF1 monoclonal antibodies to measure and visualize the concentration of secreted PRF1 via confocal microscopy. The researchers subsequently quantified the optical density of microscopy images using ImageJ 45 as a measure of secreted PRF1.
We extracted a total of 50 data points from the published study, from which 34 were used for model training and the remaining 16 were used for model validation. The experimental design in the study by Srpan et al. 26 is as follows: NK cells were stimulated under one pathway for 60 min and then washed for 15 min prior to the next round or iteration of receptor stimulation. The cells were stimulated for two rounds under one pathway and then either stimulated via the same pathway or the other pathway for the third round. We designated the data from the first two rounds of receptor stimulation for model training, and if the same pathway was stimulated again for the third round, then these data would also be assigned to the training set; otherwise, they would be assigned for model validation. For data where NK cells were stimulated for 120 min in a single well (i.e., under one pathway), we used the first half for model training and the second half for model validation. Finally, since the model predictions are units of concentration (lM) and the experimental data are signal intensity, we normalized both the model predictions and the data prior to training by calculating the fold change (FC) from the same reference time point, as done in our previous work, 19

Parameter estimation
The model parameters were estimated using a Bayesian framework, namely, the Metropolis-Hastings (MH) algorithm, 46-49 as we did previously. 19 We provide an extensive description of this approach in the published work. Briefly, the goal of this approach is to sample from the posterior probability distribution of the parameters given the data (p hjy ð Þ). Bayes' theorem provides a relationship between the posterior and prior probability densities via where p h ð Þ is our prior knowledge of the parameters, p yjh ð Þ is the data likelihood function, and p y ð Þ is the probability of the data (which is constant here since the data are given). The seven estimated model parameters are related to secretion and are briefly described in Table I. All other parameters are set to their estimated values from our previous model. 19 We used a continuous uniform prior distribution on the parameters, as there is no evidence to suggest a priori where these parameters may lie on the positive real line. Specifically, where each independent h i is an element in the arbitrary interval a i ; b i ½ . Indeed, the prior distribution is independent of h and, thus, will not play a determining role in the estimation process (see below).
The data likelihood function represents a measure of the error between the data and the model. We assume that the errors between each data point and model prediction (e i ¼ y i À M i h ð Þ) are independent and normally distributed with zero mean and identical variance r 2 , where each p y i À M i h ð Þ À Á $ N 0; r 2 ð Þ and 34 represents the number of training data. Moreover, we marginalize the noise (r 2 ) from p yjh; r 2 À Á by assuming an inverse gamma measure over r 2 (specifically, p r 2 ð Þ $ C À1 a; b ð Þ) and integrating p yjh; r 2 À Á with respect to r 2 to obtain Here, K is a constant that is independent of both y and h. Note that p yjh ð Þ achieves its maximum when y ¼ M h ð Þ, which implies that maximizing the posterior density p hjy ð Þ is proportional to minimizing the sum of squared residuals.
However, we cannot solve for p hjy ð Þ analytically due to the nonlinearities in M h ð Þ. To that end, we employ the MH algorithm to sample from the posterior distribution. Before doing so, we specify key components needed to implement the MH algorithm. 47,48 We first logtransform the parameters (i.e., calculate lnh) and estimate the parameter values in the log-space. Specifically, we use the lognormal distribution to propose a new parameter vector given the current estimate, where h Ã is the current estimate, h j ð Þ is the proposed parameter vector, j is the number of iterations of the MH algorithm, and v ¼ 0.1 is the scale parameter of the distribution. The scale parameter (v) must be carefully tuned to allow for both sufficient exploration of the (log) parameter space and maximizing the posterior density. It is important to note that the skewness and kurtosis (third and fourth moments, respectively) of the lognormal distribution grow proportionally to the exponential of the square of the scale parameter [i.e., exp v 2 ð Þ ]. As we increase the scale parameter, we exponentially increase our sampling space; however, this comes at the cost of proposing parameter vectors with low acceptance probability. Alternatively, if the scale parameter is too small, we do not sufficiently explore the parameter space and possibly may not find parameter vectors that maximize the posterior density. In our estimation, we found that when v ¼ 0.1, the MH algorithm is effective at both exploring the parameter space and maximizing the posterior density. Next, we compute the acceptance ratio (AR) at each iteration j, Since a; b > 0 is needed to satisfy the inverse gamma distribution, we set both equal to one. These two hyperparameters determine the degree of irreducible noise between the data and model predictions, where a and b affect the shape and scale of the noise, respectively. Also, note that the prior distributions do not affect the acceptance ratio since p h j ð Þ ð Þ p h Ã ð Þ ¼ 1 regardless of the parameter vectors h Ã or h j ð Þ and independent of the arbitrary intervals a i ; b i ½ for each i th component. This is a consequence of imposing a uniform prior, which implies that our knowledge of the parameters is completely determined by the data. The last term on the right-hand side is the ratio of the transition kernels, which measures the asymmetric probability of transitioning between h Ã and h j ð Þ . Given the above proposed distribution and AR, we simulated the MH algorithm for 10 000 iterations. In our estimation, the marginal posterior distribution of each parameter converges to a stationary distribution anywhere between the 2000th and 5000th iteration. Given that the MH algorithm is a stochastic optimization method, we simulated the MH algorithm 200 times using independent random initial guesses for h Ã . In addition, the uncertainty in the parameter estimates from our previous model 19 was incorporated in our current estimation of the new parameters by randomly sampling each previously fitted parameter from its corresponding distribution during every iteration of the MH algorithm. For simulations, we used the final 1000 iterations of the MH algorithm, which we take to be Indeed, inputs with a higher sensitivity index suggest that the output is sensitive to variations in the input. In our case, X i represents the kinetic parameters in our model, whereas Y 1 and Y 2 represent the amount of secreted GZMB and PRF1, respectively, after 60 min of receptor stimulation to mimic the experimental conditions from the study by Srpan et al.. 26 We drew 250 random samples for each X i (independently) using a uniform distribution on the interval 0:5E X i ð Þ; 1:5E X i ð Þ ½ , where E is the expectation operator and the distribution of each X i is given by parameter estimation (see above) or from the study by Makaryan and Finley. 19 Next, we simulated the model using these 250 samples to generate a distribution for each of the outputs Y 1 and Y 2 to then compute the total order sensitivity indices.

synNotch signaling and RNA expression model
The synthetic biology field has empowered biologists with tools for engineering novel cellular responses. In general, such molecular programming techniques provide cells with additional capabilities; for instance, Smole et al. 51 engineered novel genetic circuits in mammalian cells to respond to inflammatory signals (i.e., IL-1 b) by producing anti-inflammatory proteins. In addition, the use of synthetic biology methods with clinical applications is the subject of recent reviews. [52][53][54] The synthetic Notch (synNotch) signaling pathway, 55 in particular, can be constructed via genetic modifications to trigger a specific cell response when a specific stimulus (e.g., chemical and thermal) is present in the micro-environment. Briefly, the synNotch system includes three components: (1) the synNotch receptor, (2) a transcription factor, and (3) a plasmid vector. The extracellular domain of the synNotch receptor can be a single-chain variable fragment (scFv) designed to bind to a specific antigen, similar to chimeric antigen receptors (CARs) that target tumor-specific antigens. 16,56,57 The synNotch receptor and the transcription factor are linked together using peptide sequences that can be cleaved by constitutively expressed membrane proteases once the receptor binds to a specific ligand. 55 Then, the transcription factor becomes unchained from the cell membrane and subsequently free to bind to its promoter site and initiate gene expression.
It follows that if the synNotch receptor's extracellular domain is engineered to bind to the same ligands as CD16 and NKG2D, then this synthetic pathway will be complementary to the endogenous pathway, that is, the endogenous pathway will mediate the secretion of GZMB and PRF1, while the synthetic pathway will both increase the production of GZMB and PRF1 in addition to enhancing their secretion. To determine if this approach is indeed beneficial, we constructed a mathematical model of the synNotch receptor, based on mass action kinetics, to simulate and predict if such a system necessarily leads to the continuous secretion of GZMB and PRF1 over multiple rounds of stimulation.
We implemented an in silico synNotch system to inhibit protein activity. For inhibiting protein activity, the use of long non-coding RNAs (lncRNAs) presents new avenues for impeding protein-to-protein interactions in signaling pathways. [58][59][60][61] Recently, lncRNA pulldown assays, 61-63 when coupled with mass spectrometry, have been utilized for identifying novel RNA molecules that can bind to and sequester proteins from signaling. An advantage of proteinsequestering lncRNAs is that, once discovered, they can be sequenced and reverse-transcribed via polymerase chain reaction (PCR) to manufacture their complementary DNA strand (cDNA). 64 The cDNA can be subsequently integrated into a plasmid vector for expression under the control of the synNotch receptor. In particular, we simulate inhibition of phosphatase activity.
We also apply the synNotch system to promote production of cytolytic molecules. By incorporating a multi-cistronic plasmid downstream of the synNotch receptor, which enables the expression of two or more genes, the production of GZMB and PRF1 can be induced once the NK cell binds to a specific target ligand. This added specificity allows the researcher to designate when protein production should occur. Given that protein expression is a costly function, this strategy limits cytolytic molecule production to situations where an NK cell interacts with a target cell. Thus, the NK cell allocates its resources for cytolytic molecule production only in cases where a target cell is nearby. In contrast to constitutive over-expression techniques, this method allows the NK cell to preserve its energy for other functions.
To simplify our analysis, we assume that (1) a lncRNA binds to, and sequesters, both phosphorylated and unphosphorylated phosphatase at a rate of 1 lM Â min ð Þ À1 and (2) there are two plasmids (at fixed amounts) controlling the expression of lncRNA and the cytolytic molecules separately and both are under the control of the same transcription factor. For a detailed description of the reactions and parameters that characterize the synNotch signaling system, see Table II.
In Eq. (1), the ligand can be either Rituximab or MICA, to signal complementarily with the CD16 or NKG2D pathway, respectively. For Eq. (5), the production of RNA is nonlinear with respect to the transcription factor (TF) and follows a Hill function with plasmid affinity K M and Hill coefficient (n) equal to one. Initially, we assume no cooperativity between the TF and the plasmid; however, we relax this assumption and observe how the synNotch dynamics are affected using various values of the Hill coefficient (see the Results). We assume that the plasmid concentration (u) in Eq. (5) is at steady state prior to receptor stimulation. Rates of protein production 65 and RNA degradation 66,67 are from published data. We performed a sensitivity analysis to determine the robustness of the model predictions with respect to our assumptions on the kinetic rates.

Optimization of the synNotch system
We apply optimal control theory to determine how a synthetic pathway, the synNotch pathway, can be used to promote maximal secretion of cytolytic molecules. Such an optimization is necessary because, while the synNotch pathway may lead to an increase in the production of cytolytic molecules and the inhibition of protein activity, it certainly imposes a burden on the cell to express this system. Moreover, given that the synNotch and endogenous receptors will be competing for the same ligand, it is not immediately clear how much synNotch receptor is optimal for a given frequency of NK cell stimulation. Therefore, we considered optimizing the synthetic pathway to maximally secrete GZMB and PRF1 while using the absolute minimal amount of exogenous material. In this context, the plasmids and synNotch receptor can be considered as the controllers for the secretion of cytolytic molecules. Thus, our objective is to find the optimal values of the plasmids and synNotch receptor (controllers) such that we maximally induce secretion (performance) at a minimal cost to the cell (effort). Indeed, this is an optimization problem, which we can solve using conventional methods; 68-71 specifically, we minimized the following stochastic objective function, given the model parameters h: : Indeed, we minimize the objective function using the sample average approach (SAA). [72][73][74][75][76][77][78] Here, u; v, and R 0 represent the amount of lncRNA-coding plasmid, cytolytic molecule-coding plasmid, and the initial value of the synNotch receptor, respectively. The minimization is subject to the constraints: 0 u; v 1000 copies Â cell À1 À Á 0 R 0 10 lM ð Þ ( . We set the upper bound of the plasmid concentrations based on what is defined as a high copy number. 51,54,55 In addition, the upper bound on the initial value of synNotch receptor is based on a study 76 where Chinese hamster ovary (CHO) cells were genetically modified to maximally produce human IgG; the values ranged from 0.3-20 lM, from which we chose 10 lM. While we acknowledge the dissimilarities between the CHO and NK cells and the synNotch receptor and human IgG, it is, nevertheless, a strict upper bound that can constrain our estimation. Additionally, we solved the above objective function where the upper bound of R 0 was either 20 lM or unbound and observed no differences in the conclusion of the results. The results show that the optimality of synNotch was more sensitive to the transcriptional dynamics. The second term on the right-hand side represents the cumulative secretion of GZMB and PRF1 after 60 min of receptor stimulation, where r ¼ 1; …; N r is the number of rounds of receptor stimulation and the expectation is taken over the parameters h. Since GZMB and PRF1 are non-negative, the second term is, in fact, a maximization problem given that argmin Àf ¼ argmax f for all non-negative f . The constant q 2 0; 1 ð Þis a weight parameter that specifies how much emphasis is placed on minimizing the first term (effort) vs maximizing the second term (performance); in this study, we place an equal weight on both (i.e., q ¼ 1 2 ). To solve this optimization problem, we used the mesh-adaptivedirected search (MADS) algorithm patternsearch in MATLAB. This is a gradient-free method that attempts to locate a minimizer of the objective function by evaluating many trial points nearby the initial guess at each iteration. If some trial point near the initial guess induces a lower function evaluation, then the iteration terminates and starts again by implicitly creating a new mesh around this new incumbent point. If the algorithm cannot find a feasible point that minimizes the objective function, the mesh around the current incumbent point becomes finer and finer until a predefined threshold is reached. For our purposes, we set the mesh tolerance parameter to 10 À6 at which the algorithm terminates.

SUPPLEMENTARY MATERIAL
See the supplementary material for (File S1) computational model files (.m files with model code and .mat files with parameter values); (File S2) list of model species, reactions, and parameters (provided as .xlsx file); and (File S3) supplementary figures (provided as .pdf file).

ACKNOWLEDGMENTS
We immensely appreciate the Finley research group for critical evaluation of this manuscript as well as improvements to the model code. This work was supported a Viterbi/Graduate School Merit Fellowship (to S.Z.M).

DATA AVAILABILITY
The data that support the findings of this study are available within this article and its supplementary material.