Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Insights from a computational approach to natural killer cell activation
(USC Thesis Other)
Insights from a computational approach to natural killer cell activation
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Insights from a Computational Approach to Natural Killer Cell Activation
by
Sahak Zack Makaryan
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
Biomedical Engineering
August 2021
Copyright 2021 Sahak Zack Makaryan
ii
Dedication
This thesis is dedicated to my grandmother, who initiated my journey,
my parents and sister, who supported my journey,
and my wife, who helped me finish my journey.
iii
Acknowledgements
I want to thank my advisor Professor Stacey Finley by recalling two important
instances in my trajectory that made this work possible. During laboratory rotations, I
contacted Professor Finley to express my reluctance in joining her laboratory as I was afraid
that I would not succeed and would fail rather quickly since I did not have any prior coding
or a strong mathematical background. To my surprise, Professor Finley did not accept my
request and instead she encouraged me to join the laboratory. This was the moment I
realized that Professor Finley does not give up easily and this was exactly what I was looking
for in a mentor. This attribute of hers made all the difference for me. Thus, I joined her
laboratory immediately. Thank you, Professor Finley, for not giving up on me when it was
the easier thing to do. The second moment occurred during the end of my second year. I
walked into Professor Finley’s office for our weekly meeting to discuss my research progress.
After a few exchanges, Professor Finley told me to apply control theory to my current project.
I boldly said “Yes, absolutely! I will get it done,” without any idea of what control theory was
or how to apply it. Little did I know, she was doing me a favor. A few years later, I have come
to love the field and often approach real world problems (technical or otherwise) through
this lens. This way of thinking helped me advance in my research, lay out the trajectory of
my career and ultimately secure a position in industry months before graduation. Again,
thank you Professor Finley whether or not this was your intention.
In my personal life there are many individuals that served as a catalyst to my success.
Beginning with my nuclear family, my late father (Ovnan Makaryan) and my late paternal
grandmother (Yegisabet Etaryan) strongly encouraged me to go to college and study science.
iv
It is because of them that I started this journey, and I wish they were here to see who I have
become. They came to United States in the early 1990’s with nothing in their pockets, worked
endlessly to give everything to their children (and grandchildren) and passed on. Although I
can never repay the debt to both my grandmother and father, I hope my character,
personality and success made their trip to this side of the hemisphere worth it. I miss both
of them. My mother (Asmik Makaryan) and younger sister (Angela Makaryan) supported me
in every manner that they could, doing almost everything for me so that I do not have to and
can focus my time studying. Their sacrifice will never be forgotten. My paternal grandfather
(Sahak Makaryan) – who I am named after – alongside my uncles (Samvel Makaryan, Tatul
Shamoyan and Gagik Arakelyan) and late father taught be the importance of hard work,
resiliency and humility. Unbeknownst to me, these characteristics were precisely what I
needed to succeed in graduate school although I did not understand it at the time. My
paternal grandfather, who was fortunate enough to receive a sixth-grade education is Soviet
Armenia, told me to never belittle those less educated than myself as there is always
something to learn from everyone. My aunts (Marine Makaryan, Anahit Aslanyan and Anahit
Arakelyan) were strong influencers in my decision to continue my education; they
recognized my academic tendencies, and rightfully suggested that an undergraduate
education would leave me wanting for more – like drinking salt water. It would be very
remiss of me not to thank my younger cousins. Specifically, my cousins Tadevos Makaryan,
Harout Makaryan, Hakop Shamoyan and David Shamoyan. They believed in me even when I
was diagnosed with late-stage imposter syndrome. Here, I learned that they were counting
on me to succeed, which was an integral part of their belief that hard work pays off in this
v
world. From then on, it become clear to me that I had a role to play – whether or not I wanted
to – which was to serve as a beacon of light, helping the young ones to stay clear of the rocks.
My dear and lovely wife (Arpine Berekian), I hope I have made you proud. I am
grateful for your support through these difficult times and for putting up with my
frustrations along the way. I will never forget what you have done for me and now I will
support you as long as I am alive. The completion of this thesis is an achievement I share
alongside my wife. She did not want me to settle for anything less, realizing my potential
when I couldn’t. Moreover, I want to extend my gratitude to my in-laws (Mesrop Berekian
and Arshalouis Berekian) as well as my sister-in-law (Victoria Berekian) and her family (Alik
Piliposyan, Sofia Piliposyan and Alexa Piliposyan) for all their support and believing in me.
This work will not be what it is if it weren’t for the computational systems biology
laboratory. They were instrumental in advancing my work by telling me where I was right
and where I was wrong. Specifically, Dr. Mahua Roy, Dr. Jennifer Rohrs, Dr. Jess Wu, Ding Li,
Min Song, Colin Cess and Patrick Gelbach, were extremely helpful, kind and made the process
much more enjoyable. Also, my committee members Dr. Keyue Shen, Dr. Nicholas Graham
and Dr. Michael Khoo were absolutely necessary for providing much needed guidance to
complete this body of work. Overall, if I had to do it all over again, I would not change a thing.
vi
TABLE OF CONTENTS
Dedication .................................................................................................................................................................... ii
Acknowledgements ................................................................................................................................................ iii
List of Tables .............................................................................................................................................................. ix
List of Figures ............................................................................................................................................................. x
Abstract ........................................................................................................................................................................ xi
Chapter 1: Introduction ......................................................................................................................................... 1
1.1 Natural killer (NK) cells ................................................................................................................... 1
1.2 The signaling network of NK cells ............................................................................................... 2
1.2.1 Dynamics of membrane proximal molecules ....................................................... 3
1.2.2 Dynamics of membrane distal molecules .............................................................. 4
1.3 Systems biology ................................................................................................................................... 6
1.4 Dissertation outline ........................................................................................................................... 7
Chapter 2: Enhancing network activation in natural killer cells: predictions from in silico
modeling ............................................................................................................................................... 9
2.1 Abstract ................................................................................................................................................... 9
2.2 Introduction ....................................................................................................................................... 10
2.3 Methods ............................................................................................................................................... 14
2.3.1 Model Construction ...................................................................................................... 14
2.3.2 Data collection and processing ................................................................................ 16
2.3.3 Parameter estimation .................................................................................................. 18
2.3.4 Construction of magnitude of network activation metric ........................... 22
2.3.5 Clustering and principal component analysis ................................................... 24
2.4 Results .................................................................................................................................................. 24
2.4.1 Model of NK cell signaling matches experimental data ................................ 24
2.4.2 Baseline network activation is greatest when CD16 is stimulated ......... 29
2.4.3 Receptor characteristics significantly influence network activation ..... 30
vii
2.4.4 Co-stimulation of CD16 and NKG2D potently activates the network ..... 35
2.4.5 In silico perturbation highlight the role of phospho-receptors and
phosphatases in enhancing network activation .......................................................... 37
2.5 Discussion ........................................................................................................................................... 40
2.6 Acknowledgements ......................................................................................................................... 44
2.7 Appendix ............................................................................................................................................. 44
Chapter 3: An optimal control approach for enhancing natural killer cells’ secretion of
cytolytic molecules ........................................................................................................................ 58
3.1 Abstract ................................................................................................................................................ 58
3.2 Introduction ....................................................................................................................................... 59
3.3 Methods ............................................................................................................................................... 63
3.3.1 Construction of NK cell degranulation model ................................................... 64
3.3.2 Data collection and processing ................................................................................ 66
3.3.3 Parameter estimation .................................................................................................. 68
3.3.4 Information-theoretic sensitivity analysis ......................................................... 72
3.3.5 synNotch signaling and RNA expression model .............................................. 73
3.3.6 Optimization of synNotch system .......................................................................... 76
3.4 Results .................................................................................................................................................. 79
3.4.1 NK cell degranulation model can reproduce experimental observations
.......................................................................................................................................................... 79
3.4.2 Inhibition of pSHP maximizes PRF1 and GZMB secretion in silico ......... 82
3.4.3 The optimal synNotch system depends on the number of rounds of
stimulation ....................................................................................................................... 85
3.4.4 The optimal CD16- and NKG2D-synNotch pairs show qualitative
differences ........................................................................................................................ 90
3.4.5 Transcription factor cooperativity impacts the performance of
synNotch ............................................................................................................................ 95
3.4.6 The synthetic-endogenous pathway pairs promote a biphasic response
to stimulus ........................................................................................................................ 99
3.5 Discussion ........................................................................................................................................ 102
3.6 Acknowledgements ...................................................................................................................... 108
3.7 Appendix .......................................................................................................................................... 108
Chapter 4: A population-balance model of natural killer cell and tumor cell interactions 125
viii
4.1 Abstract ............................................................................................................................................. 125
4.2 Introduction .................................................................................................................................... 126
4.3 Methods ............................................................................................................................................ 129
4.3.1 Population-balance modeling ............................................................................... 129
4.3.2 NK cell population ...................................................................................................... 131
4.3.3 Tumor cell population .............................................................................................. 134
4.3.4 Complex population .................................................................................................. 136
4.4 Results ............................................................................................................................................... 139
4.4.1 The initial rate of tumor growth is slower if the tumor selects for
evasion ............................................................................................................................ 139
4.5 Discussion ........................................................................................................................................ 141
4.6 Acknowledgements ...................................................................................................................... 143
Chapter 5: Conclusions ..................................................................................................................................... 144
5.1 Summary .......................................................................................................................................... 144
5.2 Future directions .......................................................................................................................... 146
5.3 Conclusion ....................................................................................................................................... 148
References: ............................................................................................................................................................ 149
ix
List of Tables
Table 4-1. List of cellular parameters.........................................................................................................138
Table 4-2. List of molecular parameters....................................................................................................138
x
List of Figures
Figure 1-1. NK cell signaling model...................................................................................................................5
Figure 2-1. Model schematic..............................................................................................................................15
Figure 2-2. Model calibration and validation.............................................................................................27
Figure 2-3. Baseline network activation of individual receptors......................................................29
Figure 2-4. Network activation as a function of ligand stimulation.................................................31
Figure 2-5. Receptor co-stimulation..............................................................................................................35
Figure 2-6. Perturbations to the stimulatory network..........................................................................39
Figure 3-1. Natural killer cell signaling model schematic.....................................................................65
Figure 3-2. Model training and validation...................................................................................................81
Figure 3-3. Sensitivity analysis shows that inhibition of pSHP activation is most influential
for GZMB and PRF1 secretion.....................................................................................................84
Figure 3-4. Optimal strategies depend on the number of rounds of stimulation.......................87
Figure 3-5. Transcription factor cooperativity impacts CD16-synNotch performance...........97
Figure 3-6. Long-term behavior of the combined synthetic-endogenous pathway...............100
Figure 4-1. Model schematic of NK- and tumor-cell interactions...................................................130
Figure 4-2. Initial distribution of NK cells.................................................................................................139
Figure 4-3. Probability measure on Ω
!
......................................................................................................140
Figure 4-4. Initial dynamics of NK and tumor cells...............................................................................141
xi
Abstract
The interactions between natural killer (NK) cells – a subset of effector immune cells
– and tumor cells has gained traction over the past decade due to the innate tumor cell killing
ability of NK cells. Indeed, researchers have developed many therapeutic strategies, ranging
from antibody-directed methods (e.g., bi-specific engagers) to genetic modifications that
promote NK cell activation upon contact with cancer cells. Unlike conventional
chemotherapy, the efficacy of such cell-mediated therapies relies heavily on the NK cell’s
ability to eliminate cancer cells in an effective and successful manner. There are a couple of
non-negligible drawbacks with this approach: (1) NK cells are heterogenous both in their
killing potential and in their expression of stimulatory receptors that mediate cancer cell
killing and (2) the activation of NK cells through its stimulatory signaling network is
nonlinear with respect to its input. Moreover, finding optimal therapeutic strategies requires
a quantitative understanding of the aforementioned issues, and without such a framework
our therapies will almost surely be sub-optimal. In this work, I address the above concerns
by constructing and simulating a set of mathematical models that describes the inter- and
intra-cellular signaling dynamics between NK cells and tumor cells. The models are applied
to predict (1) which stimulatory receptors – alone or in combination – most potently activate
NK cell signaling network, (2) which strategies enhance NK cell activation with respect to
the availability of input, (3) which intracellular signaling species most strongly controls and
regulates NK cell activation, (4) how to perturb the signaling network to unbridle cell
activation, (5) how to modify NK cells – via synthetic biology approaches – such that they
secrete more cytolytic molecules upon receptor stimulation and (6) whether the density of
xii
stimulatory receptors on the cell surface of NK cells is strong determinant of tumor growth
inhibition. Ultimately, the application of such mathematical models can inform both clinically
relevant therapeutic strategies as well as a deeper understanding of NK cell signaling and
activation.
1
Chapter 1: Introduction
1.1 Natural killer (NK) cells
Natural killer (NK) cells are immune cells that belong to the host’s innate immune
system
1,2
. As their name suggests, NK cells are capable of killing diseased cells – including
both cancer cells as well as virally-infected cells – which makes them an attractive platform
for immune cell-mediated therapies (i.e., immunotherapies). NK cells possess two primary,
but not mutually exclusive, mechanisms mediating target cell lysis
3
: (1) the activation of key
stimulatory receptors promoted by agonistic antibodies or stimulatory antigens on the
target cell surface and (2) the blocking of inhibitory receptors facilitated by antagonistic
antibodies (e.g., checkpoint blockades) as well as an absence of inhibitory antigens on the
surface of target cells. In other words, NK cells can be activated directly via stimulatory
ligands and indirectly via the absence of inhibitory ligands. Indeed, it is highly likely that both
mechanisms of target cell killing occur in vivo. Irrespective to the mechanism of action, the
activation of NK cells triggers the release of perforin-1 (PRF1) and granzyme B (GZMB),
which mediates the formation of pores on the target cell surface allowing for the GZMB to
infiltrate the intracellular compartment of the target cell and facilitate apoptosis (i.e.,
programmed cell death). This intrinsic ability of cell killing has prompted many researchers
to investigate how NK cells can be directed to eliminate specific disease-promoting cells that
drive a particular pathology.
The clinical applications of NK cell-based therapies have gained much attention from
both scientists and physicians alike, especially in the interdisciplinary field of cancer
2
immunotherapy. In fact, patients suffering from cancers with a hematological origin have
benefited greatly from such therapies
4,5
. Unfortunately, the same level of success has not
been seen for patients with solid tumors. This disparity in patient outcomes is thought to be
due to the dissimilarities in the tumor microenvironment between such types of pathologies.
Research has shown several mechanisms – both individually and collectively – that may
explain the suboptimal performance of NK cell-based immunotherapies, ranging from the
mechanical properties of the tumor
6–9
, the chemical composition of the extracellular
environment
10–15
and the functional impairment of NK cells
16–20
. This work is focused on the
latter of the aforementioned processes.
1.2 The signaling network of NK cells
NK cells distinguish between normal and diseased cells through receptors on their
cell surface
21,22
. Through a diverse repertoire of stimulatory and inhibitory receptors, NK
cells are able to recognize if a target cell needs to be eliminated or not. In particular, normal
cells express major histocompatibility molecules of class I (MHC-I) on their surface, which
binds to NK cell inhibitory receptors (e.g., CD94/NKG2A). Indeed, such inhibitory antigens
are necessary to prevent NK cells from attacking the host. Diseased cells can lose their ability
to sufficiently express MHC-I, effectively rendering them as potential targets since they may
not be able to preclude NK cell activation upon contact. In addition, transformed cells can
upregulate the expression of stress-induced antigens (e.g., MHC-I polypeptide-related
sequence A (MICA)), which functionally serve as a signal to NK cells to promote cell
activation and target cell lysis. The presence of antibodies on the surface of target cells can
3
direct NK cells to kill by activating stimulatory receptors that bind to the constant region of
such antibodies. Collectively, these mechanisms enable NK cells to designate the contacted
cell as either a friend or foe.
1.2.1 Dynamics of membrane proximal molecules
The structural nature of surface receptors largely determines whether or not the NK
cell will ultimately activate and secrete cytolytic molecules (e.g., perforin-1 (PRF1),
granzyme B (GZMB)). In particular, the surface receptors can be covalently bound to
intracellular adaptor molecules located near the membrane. These adaptor molecules are
indispensable mediators of cell activation. For instance, CD16 is an Fc𝛾 receptor that binds
to the Fc portion of antibodies. Moreover, CD16 is co-expressed with either the CD3𝜁 or 𝛾-
chain molecules
23
. These species contain unique amino acid sequences, known as
immunoreceptor tyrosine-based activation motifs (ITAMs), that attract a family of kinases
when they become phosphorylated. In turn, these kinases then begin to phosphorylate many
membrane-bound scaffolding proteins such as LAT to bring together the signalosome – a
complex of proteins that increases the local concentration and signaling activity of its
components. This positive feedback mechanism promotes the self-assembly and clustering
of membrane-proximal signaling molecules near the receptor and thus initiates signal
transduction
24,25
. The absence of ITAMs in CD3𝜁 nullifies the stimulatory nature of the CD16
receptor, implying these sequences of amino acids are required for cell activation. The ITAM
sequences, however, are not the only sequences capable of generating such self-
organization. The inhibitory receptor CD94/NKG2A contains immunoreceptor tyrosine-
4
based inhibitory motifs (ITIMs) that attract phosphatases (e.g., SHP-1, SHP-2) to
dephosphorylate and inhibit the activation of the signalosome. In general, such motifs
determine whether the receptor pathway is stimulatory or inhibitory in nature.
1.2.2 Dynamics of membrane distal molecules
In addition to the receptor or adaptor molecule determining the qualitative nature of
the pathway, the cytosolic signaling species influence the magnitude of cell activation or
inhibition. Once the signalosome is assembled, the Src family of kinases (SFK) phosphorylate
and mediate the activation of key signaling molecules (Figure 1-1). Specifically, SFK can
activate members of the Syk family of kinases (e.g., ZAP70) in addition to the scaffolding
protein LAT. The phosphorylation of LAT recruits the proteins PLC𝛾, Vav and SLP76 to the
membrane where they can become activated by SFK. The molecule PLC𝛾, once activated,
ultimately facilitates the release of calcium ions from the endoplasmic reticulum, which
polarizes the cell membrane, thereby enabling the docking and exocytosis of secretory
vesicles
26–28
. The activation of SLP76 and Vav induces actin remodeling such that the cell
polarizes its actin filaments towards the local cluster of activated receptors. This directs the
secretory vesicles to the immunological synapse between the NK cell and the target cell.
Moreover, LAT can attract the kinase PI3K, which – upon activation via SFK – leads to the
activation of Akt, an important signaling molecule for cell survival. Similarly, LAT can recruit
and promote the activation of the Ras-Raf-MEK-Erk pathway, which is correlated with cell
proliferation. Interestingly, SFK can mediate its own inhibition indirectly by activating the
phosphatase SHP in an incoherent feed-forward mechanism. While counterintuitive, this
5
signaling motif is not uncommon in biological systems
29
; it helps to prevent over-activation.
Overall, the dynamic interplay between SFK and SHP impacts the extent of species activation
and thus cell activation.
Figure 1-1 | NK cell signaling model. Stimulatory receptors (green) activate the signaling
network via phosphorylation of SFKs (pSFK) whereas inhibitory receptors (red) inhibit cell
activation via activation of the phosphatase SHP. Black arrows, stimulation; red crossbars,
inhibition.
6
1.3 Systems biology
The dynamics of cell signaling pathways in general are difficult to disentangle by
experimentation alone, since the time evolution of the molecular species are almost always
nonlinear, coupled to one another with positive or negative feedback loops (or sometimes
both). Excitingly, systems biology – a subfield that utilizes mathematical modeling to
describe biological processes – can be leveraged to obtain a more quantitative understanding
of such networks. Moreover, if experimental data or prior knowledge are available, the
model predictions can be constrained using conventional model fitting approaches to match
experimental observations. These so-called data-driven models are especially useful in
clinical and translational studies as they yield predictions that are physiologically realistic.
Indeed, such a feature that can be helpful to clinicians when deciding on, say, chemotherapy
dosing strategies since the physician will not have all the time and resources to test all
possibilities. A mathematical model, on the other hand, can be simulated thousands of times,
over an exhaustive range of parameters and initial conditions, in a fraction of the time to do
a single experiment. Still, realistic predictions notwithstanding, mathematical models of
biological systems force us to explicitly delineate our assumptions about the underlying
system, and therefore generate hypotheses to test and falsify these assumptions – which can
help expand our knowledge, or re-adjust our prior knowledge, of the system being studied.
In total, our understanding of biology benefits the most when there is a union between
experimentation and mathematical modeling.
7
1.4 Dissertation outline
In this work, I present a series of mathematical models that predict the signaling
dynamics of NK cells, and apply these models to inform strategies that enhance NK cell
activation in the context of cancer cell killing.
In Chapter 2, I construct an experimentally validated, mechanistic model of three
natural killer cell signaling pathways (CD16, 2B4 and NKG2D) to provide insight on how to
enhance activation of key cytotoxicity-mediating species. Here, I show that the receptor
concentrations and their rate of activation influence the maximum level of species activation
and the amount of stimulus needed to attain maximal activation, respectively. Furthermore,
it will be demonstrated that the co-stimulation of CD16 and NKG2D most effectively activates
the signaling species compared to other combinations. Lastly, I perturb the model to
highlight how the availability of ligand determines which strategy is most potent in
activating the signaling species.
In Chapter 3, I extend the model described in Chapter 2 to include the dynamics of the
cytolytic molecules GZMB and PRF1. Moreover, I fit the model predictions to published
experimental data by estimating key kinetic rate constants using a Bayesian perspective.
Furthermore, I employ an information-theoretic sensitivity analysis to demonstrate which
nodes in the CD16 and NKG2D signaling networks greatly influence the amount of secreted
GZMB and PRF1. With this knowledge at hand, we expand the model further to include a
synthetic pathway to (1) promote the production of GZMB and PRF1 and (2) to increase the
amount the secreted GZMB and PRF1 by disinhibiting signal transduction. Finally, I apply
8
optimal control theory to the synthetic pathway to identify the minimal amount of
exogenous material that maximizes the amount of secreted GZMB and PRF1.
In Chapter 4, I construct a mathematical model that details how NK cells interact with
tumor cells and how the cell types change their receptor concentrations as a result of such
interactions. In effect, I use a system of partial differential equations (PDEs) to model the
dynamics. Furthermore, I approximate the solution of the system of PDEs using the finite-
volume method, and I demonstrate how the initial receptor distribution of NK cells and
which properties of the tumor cells accelerate tumor growth. Altogether, my work focuses
on using mechanistic modeling to predict the dynamics of NK cell signaling and inform
strategies that maximize their effectiveness in the context of tumor cell killing.
9
Chapter 2: Enhancing network activation in natural killer cells:
predictions from in silico modeling
Portions of this chapter are adapted from:
Sahak Z. Makaryan and Stacey D. Finley.
Integrative Biology (2020), Vol. 00, No. 00: 1 – 13
2.1 Abstract
Natural killer (NK) cells are part of the innate immune system and are capable of
killing diseased cells. As a result, NK cells are being used for adoptive cell therapies for cancer
patients. The activation of NK cell stimulatory receptors leads to a cascade of intracellular
phosphorylation reactions, which activates key signaling species that facilitate the secretion
of cytolytic molecules required for cell killing. Strategies that maximize the activation of such
intracellular species can increase the likelihood of NK cell killing upon contact with a cancer
cell and thereby improve efficacy of NK cell-based therapies. However, due to the complexity
of intracellular signaling, it is difficult to deduce a priori which strategies can enhance species
activation. Therefore, we constructed a mechanistic model of the CD16, 2B4 and NKG2D
signaling pathways in NK cells to simulate strategies that enhance signaling. The model
predictions were fit to published data and validated with a separate dataset. Model
simulations demonstrate strong network activation when the CD16 pathway is stimulated.
The magnitude of species activation is most sensitive to the receptor’s initial concentration
and the rate at which the receptor is activated. Co-stimulation of CD16 and NKG2D in silico
10
required fewer ligands to achieve half-maximal activation than other combinations,
suggesting co-stimulating these pathways is most effective in activating the species. We
applied the model to predict the effects of perturbing the signaling network and found two
strategies that can potently enhance network activation. When the availability of ligands is
low, it is more influential to engineer NK cell receptors that are resistant to proteolytic
cleavage. In contrast, for high ligand concentrations, inhibiting phosphatase activity leads to
sustained species activation. The work presented here establishes a framework for
understanding the nonlinear aspects of NK cell signaling and activation.
2.2 Introduction
Natural killer (NK) cells are immune cells that can eliminate cancer cells upon cell
contact
1,2,30
. NK cells express a repertoire of stimulatory receptors that mediate the release
of cytotoxic chemicals when stimulated by antibodies or by cells that express stimulatory
ligands. The activation of such receptors induces intracellular signaling through a cascade of
phosphorylation reactions, which ultimately leads to NK cell activation, the secretion of
cytolytic molecules and cancer cell death. This innate ability for cancer cell elimination has
spurred an interest in research
30–32
to better understand NK cell activation. In vitro
studies
26,33–35
suggest a strong correlation between the activation of key signaling species
and NK cell activation via target cell killing assays such as
51
Cr-release assays as well as
cytokine production via ELISA. Thus, we hypothesize that enhancing the activation of key
signaling species can proportionally enhance cancer cell killing and thereby improve patient
outcomes in the clinic. Given that cancer cell killing is initiated via activation of NK cell
11
stimulatory receptors, it is important to understand how the signal propagates and activates
the downstream species that contribute to target cell killing. Therefore,
researchers
1,2,30,31,34,36
have studied NK cell signaling and reported which species are
activated downstream of the stimulatory receptors. Such findings are crucial in
understanding how NK cell activation proceeds on the molecular level.
However, due to the natural complexity and nonlinearity underpinning intracellular
signaling, it is difficult to deduce how NK cell signaling can be modulated to enhance species
activation. Mathematical models are valuable in these contexts in that they enable us to
untangle such complicated system behavior and predict the system’s response to a wide
variety of perturbations
37–43
. For example, work by Das
38
demonstrated how receptor-ligand
interactions impact NK cell activation and the various NK cell responses induced by strong
and weak stimulatory ligands. Mesecke and colleagues
39
showed that the physical
association of Src family kinases (SFK) with a stimulatory receptor is essential for NK cells
to promote a cytotoxic response and that the activation of the signaling species Vav
correlates with P815 tumor cell killing. Nevertheless, the question of which strategies
enhance NK cell signaling (and why) remains open. Additionally, the previous models did
not determine which molecular perturbations or which pathways should be co-stimulated
to optimally activate the NK stimulatory network.
Here, we developed a molecularly detailed, experimentally validated mechanistic
model of NK cell signaling to address the above questions. The CD16, 2B4 and NKG2D
stimulatory pathways were modeled in this study as these pathways contribute to target cell
lysis in different ways
31,44,45
. CD16 is an Fc receptor that binds to the constant region of
antibodies. This implicates CD16’s activation in antibody-dependent cell-mediated
12
cytotoxicity. Its cytoplasmic domain is associated with CD3ζ, which contains
immunoreceptor tyrosine-based activation motifs (ITAM). 2B4 is part of the signaling
lymphocytic activation molecule (SLAM) family of receptors, and its cytoplasmic tail
contains four immunoreceptor tyrosine-based switch motifs. The ligand for 2B4, CD48, is
expressed by cells of hematopoietic origin. This suggests 2B4 may play a role in regulating
hematopoietic processes. NKG2D belongs to the family of C-type lectin-like receptors. It
associates with the adaptor protein DAP10, which has an activation motif that is similar to
the CD28 T cell co-receptor. NKG2D binds to ligands typically expressed by cells that have
undergone transformation, which implicates this receptor in the elimination of tumors. We
focus on the CD16, 2B4 and NKG2D pathways as these receptor pathways are well studied
with respect to NK cell signaling and involve different intracellular signaling components.
Other natural cytotoxicity receptors (NCRs) (e.g., NKp30, NKp46, NKp44) could have been
included in this study; however, many of these receptors share similar intracellular signaling
domains as CD16. Specifically, NKp30 and NKp46 signal via ITAM-containing proteins such
as CD3ζ, and NKp44 signals through DAP12, which also contains a single ITAM motif
3
. Thus,
to avoid redundancy in modeling similar intracellular signaling pathways, we excluded the
NCRs from this work. This allows us to study NK cell signaling via three pathways that are
not only stimulated under different physiological scenarios but also signal via different
intracellular proteins containing different activation motifs.
Ligand binding to the CD16, 2B4 and NKG2D receptors initiates intracellular
signaling. The PI3K-Akt, SLP76-Vav-Erk and PLCγ networks are all activated upon
stimulation of CD16, 2B4 and NKG2D
3
. In NK cell biology, activation of the PI3K-Akt signaling
axis is correlated with cell survival, while Erk activation is correlated with cell
13
proliferation
26,34,35
. SLP76 and Vav activation are necessary for actin remodeling and the
formation of the immunological synapse
46
. In addition, the knockdown of SLP76 and Vav via
siRNA in NK cells has been correlated with a dampened release of intracellular calcium ions
upon stimulation of NKG2D and 2B4 as well as significantly lower percent of specific lysis
against P815 tumor cells
26,34,46
. The knockdown of these molecules is also correlated with
lower production and secretion of the cytokines IFN-γ and MIP-1α. PLCγ activation induces
the release of intracellular calcium ions, which subsequently contributes to cell activation,
and correlates with cytokine production and target cell killing via
51
Cr-release assay
32
.
Indeed, activation of these signaling intermediates is necessary to activate and induce a
cytotoxic response in NK cells. By studying the initial dynamics of the signaling species, we
can predict NK cell responses to receptor stimulation. Although most of the downstream
reactions are common between the pathways, there are some subtle differences. For
example, 2B4 does not induce Akt phosphorylation
34,36
. Additionally, 2B4 and NKG2D
specifically lead to phosphorylation of the Y113 and Y128 sites on SLP76, respectively, while
CD16 induces phosphorylation of both sites
34
. Also, CD16 induces ZAP70 and LAT activation,
while 2B4 and NKG2D do not. Thus, these pathways are interconnected, and understanding
the dynamics of the concentrations of the molecular species involved in the signaling
pathways requires in-depth analyses.
In the present study, we use mathematical modeling to characterize and compare the
signaling dynamics of the CD16, 2B4 and NKG2D pathways with respect to their magnitude
of activation of the network. Furthermore, we identify which signaling species and
parameters influence the magnitude of network activation and which combinations of
receptor co-stimulation most potently activate the network. In silico perturbations of the
14
stimulatory network demonstrate the strategies that effectively increase network activation,
including which signaling species to target and how to modify the species. In total, the model
predictions can be used for engineering NK cells with enhanced signaling, which may
improve target cell killing.
2.3 Methods
2.3.1 Model Construction
We constructed an ordinary differential equation (ODE) model to predict the
concentrations of the molecular species in the CD16, 2B4 and NKG2D pathways in NK cells
(Figure 2-1). The list of model species, reactions and parameters are provided in
Supplementary File S2-1. The rates of the biochemical signaling reactions were represented
using Michaelis–Menten reactions. Arriving at the current model structure was an iterative
process where we fitted several model types (e.g., Michaelis–Menten kinetics versus mass
action kinetics, including versus excluding the phosphatases, including versus excluding
nonspecific decay rate) to the experimental data and selected the model structure that
generated the lowest error between the model predictions and experimental data. We
constructed the model using BioNetGen
47
and simulated it in MATLAB (MathWorks).
The final model contains 83 parameters and 36 species, including the three NK cell
receptors. Each receptor binds to its ligand and forms a receptor-ligand complex that allows
the receptor to become phosphorylated by basally active SFK. Then, the ligand-bound
phosphorylated receptor serves as the catalyst for converting SFK from a basally active state
15
to a fully active state (pSFK). Fully active SFK mediates the phosphorylation (activation) of
a number of downstream signaling species, including LAT, ZAP70, PLCγ, Vav, SLP76, Akt, and
the phosphatases SHP and SHIP. Moreover, the stimulation of 2B4 can lead to activation of
the phosphatases independent of pSFK
3
. Phosphorylated ZAP70 promotes activation of LAT.
The inhibitory species, phosphatases SHP and SHIP, provide negative feedback to prevent
overactivation
48
. The catalysts for Erk phosphorylation are the phosphorylated forms of
SLP76 and Vav. These species are upstream inputs to the MAPK pathway
38,39
.
The initial concentrations of the species in the system were extracted from the
literature
49–52
, where their values were measured in resting primary NK cells using
Figure 2-1 | Model schematic. Reaction network for three stimulatory receptors expressed on
the surface of NK cells: CD16, 2B4 and NKG2D. These receptors promote signaling species that
mediate NK cell activation: SFK, Erk, Akt and PLCγ. Arrows indicate stimulation, while red
crossbars indicate inhibition. Orange arrows are specific to the NKG2D pathway; blue, CD16
pathway;purple,2B4pathway;black,allpathways.
16
quantitative mass spectroscopy. These values are assumed to be at steady state since the NK
cells were preincubated in cell culture media containing no stimulatory ligands for CD16,
2B4 or NKG2D. Given the timescale of our analysis (60 minutes following receptor
stimulation), we consider the synthesis of the signaling molecules to be negligible. In fact,
the expression of CD16 and NKG2D on the surface of primary NK cells, as measured by flow
cytometry, decreases after receptor stimulation and remains unchanged in the absence of
stimulatory ligands
53
. Also, there is no data to suggest the signaling molecules significantly
change in concentration within a 60-minute time interval in the absence of a stimulus. In
total, we consider the phosphorylation and dephosphorylation reactions of the signaling
species as well as a nonspecific decay term for the phosphorylated species in the system to
account for degradation, dilution and other mechanisms of disappearance
26,32–34,54
.
2.3.2 Data collection and processing
We trained the mathematical model using experimental data extracted from the
literature. The raw data and our data processing procedure are provided in Supplementary
File S3. To control for variations in the experimental conditions, we only used data from
published studies where (1) the antibodies used for CD16, 2B4 and NKG2D stimulation were
of the same concentration (10 𝜇g/ml) and from the same vendor, and (2) the cell types used
in these studies were primary NK cells. Immunoblot images from these published studies
were analyzed and processed using ImageJ
55
. Specifically, ImageJ provides a measure of the
optical density for any predefined rectangular space of an image in grayscale, where the
estimated optical densities range from 0 to 225 (black to white, respectively). Protein bands
17
in western blots were analyzed to estimate their optical density. To control for immunoblot
variations across the experiments, we subtracted the optical density measurement of the
western blot gel background from the optical density measurements of all protein bands in
the same gel. Furthermore, for a single protein, the optical density measurement of the
zeroth time point was also subtracted from the optical density measurements of the
remaining time points, making the initial time point a true zero. This procedure, which
follows the documented ImageJ usage protocol, standardizes the experiments for
comparison and controls for the background and zeroth time point measurements. In total,
the model was trained to 64 data points. Additionally, the model was validated against 32
data points. The signal intensity rQ
"#
!
t of a given phosphorylated species (pX) at the 𝑗
$%
time
point is calculated as:
Q
"#
!
=
OD
"#
!
OD
&'($)'*
!
where OD
"#
!
and OD
&'($)'*
!
are the optical density values of the phosphorylated species and
a loading control, respectively, at the 𝑗
$%
time point. Furthermore, the signal intensity was
normalized to a single (reference) time point xQ
"#
"#$
y by calculating the percent change
x%∆
"#
y:
%∆
"#
!
=
Q
"#
!
−Q
"#
"#$
Q
"#
"#$
×100%.
18
With normalization, the signal intensity is initially zero, so the percent change at 𝑡 =
0 is −100%. For each time-course dataset, we refrained from normalizing the experimental
data to the time point where the maximum optical density was observed. This would
introduce bias in our data processing procedure as we would be assuming the true maximum
of species phosphorylation occurred at the same time point the data were collected.
Similarly, we did not normalize the data using the initial time point measurement because
the signal intensity measurement there is zero, and we cannot divide by zero. Instead, we
normalized each time-course dataset to the time point where the median optical density was
observed (see Supplementary File S3 for details) to allow the model to predict when the true
maximum of species phosphorylation may occur. Thus, the reference (normalization) time
point almost always corresponded to the 10-minute time point. To clarify, we normalize both
the experimental data and the model simulations at the same time points to effectively
compare their values during parameter estimation. Also, since we are analyzing a single
immunoblot image for the phosphorylation of an individual species, the data contain a single
replicate for each time point. However, the authors of those experimental studies stated the
immunoblot images were representative of three independent experiments, and we treat the
single replicate as the expected outcome.
2.3.3 Parameter estimation
We assume the parameters are random variables and consequently we estimate their
distribution using a Bayesian perspective
56
, where we maximized the posterior density
x𝑝(𝜃|𝑦)y the parameters (𝜃) given the data (𝑦) via the Metropolis–Hastings algorithm
56,57
.
19
In particular, we are estimating 83 parameters using 64 data points (training dataset). While
it is possible to estimate all 83 parameters using exactly 83 data points (or even more, up to
the total number of data points available), we could possibly overfit the data. That is, we may
produce a model that has low bias (small training error) but high variance (large validation
error). The resulting model would have low predictive power due to its high variance. We
balance the bias/variance tradeoff by using a slightly larger validation dataset at the expense
of a slightly smaller training dataset.
Briefly, Bayes’ theorem describes the relationship between the posterior and prior
probabilities via the formula:
𝑝(𝜃|𝑦) =
𝑝(𝑦|𝜃)𝑝(𝜃)
𝑝(𝑦)
∝𝑝(𝑦|𝜃)𝑝(𝜃),
where 𝑝(𝑦|𝜃) represents the data likelihood function, 𝑝(𝜃) is our prior knowledge on 𝜃 and
𝑝(𝑦) is the probability of the data (which is constant here since the data is given). The
parameters in the present model are 𝑘
+,$
and 𝐾
-
, constants that represent the catalytic rate
and the half-saturating constant in enzymatic reactions, respectively. The probability
distribution of these two types of parameters has been observed to follow a lognormal
distribution and are extensively overlapping for a wide range of eukaryotic enzymes
52
. The
𝐾
-
parameter can be decomposed into
.
&'(
/.
)$$
.
)*
; however, this representation would increase
the parameter search-space and over-parameterize the model. Moreover, without complete
information of the distribution of each 𝑘
'(
and 𝑘
'00
rate constant for each reaction in the
model, it is unclear where their values would lie on the positive real line. To avoid
20
complexity, we assumed the prior probability of the model parameters are independent and
drawn from an identical lognormal distribution (IID) with 𝜇
1
= 1 and 𝜎
1
= 2. That is:
𝑝(𝜃)=𝑝(𝜃
2
,…,𝜃
34
) =𝑝(𝜃
5
)
34
562
.
The data likelihood function captures the error, 𝜖, between the model predictions and
the experimental data via 𝜖 =𝑦−ℳ(𝜃). We assume the error follows a normal distribution
centered at zero with some noise (𝜎
7
). Specifically,
𝑝(𝑦|𝜃,𝜎
7
) ≝𝑝(𝜖)=𝑝x𝑦−ℳ(𝜃)y ~ 𝑁(0,𝜎
7
),
where ℳ(𝜃) is the model prediction. In addition, we marginalized out the noise from
𝑝(𝑦|𝜃,𝜎
7
) by assuming an inverse gamma distribution over 𝜎
7
(with hyperparameters 𝛼 =
2 and 𝛽 =1) and integrating 𝑝(𝑦|𝜃,𝜎
7
) with respect to 𝜎
7
to attain:
𝑝(𝑦|𝜃)= 𝑝(𝑦|𝜃,𝜎
7
)
8
1
𝑝(𝜎
7
) d𝜎
7
.
We note that the density of the data likelihood 𝑝(𝑦|𝜃) is at its maximum when 𝑦 =ℳ(𝜃) is
centered at zero. Therefore, maximizing the posterior density is equivalent to minimizing
the error between the model prediction and the experimental data. The prior 𝑝(𝜃) serves as
a penalty term that penalizes parameter values that are at the tail ends of the distribution
21
and rewards values closer to the mode of the distribution. We cannot solve for 𝑝(𝜃|𝑦)
analytically since ℳ(𝜃) in the data likelihood is a nonlinear operator. Instead, we employ
the MH algorithm to sample from the posterior distribution, which is the target distribution.
First, we initialize 𝜃
∗
by randomly sampling from 𝑝(𝜃). Then, we sample another vector 𝜃
(5)
from a proposal distribution 𝑝x𝜃
(5)
𝜃
∗
,𝜈
7
y ~ Lognormal(ln(𝜃
∗
),𝜈
7
), where 𝑖 is the number
of iterations of the MH algorithm and we fixed 𝜈
7
= 0.1. We chose the lognormal distribution
for the proposal since it has proper support over 𝜃. Then, we compute the acceptance ratio
(AR) for each iteration, where:
AR=
𝑝x𝑦𝜃
(5)
y
𝑝(𝑦|𝜃
∗
)
𝑝x𝜃
(5)
y
𝑝(𝜃
∗
)
𝑝x𝜃
∗
𝜃
(5)
,𝜈
7
y
𝑝x𝜃
(5)
𝜃
∗
,𝜈
7
y
.
Also, we randomly sample a number λ ∼ Unif(0, 1) and compare the acceptance ratio
to λ for each iteration 𝑖. If the error between the model prediction and the experimental data
is smaller for 𝜃
(5)
than 𝜃
∗
, then 𝑝x𝑦𝜃
(5)
y>𝑝(𝑦|𝜃
∗
) and AR > λ almost surely. Thus, we
accept the parameter vector 𝜃
(5)
by setting 𝜃
∗
=𝜃
(5)
and continue this process for a fixed
number of iterations. The prior distribution remains fixed over all iterations while the
proposal distribution re-centers around parameters 𝜃
∗
that minimize the error between the
model and the data.
Based on our simulations, the algorithm almost always converges to a minimum after
∼ 5000
th
iteration. Thus, we simulated the algorithm for 10,000 iterations to ensure the
algorithm converges before termination. The first several thousand iterations of the MH
algorithm serve to maximize the posterior density, thereby converging our initial estimate
22
of the posterior distribution closer to the true posterior distribution. This is known as the
burning-in phase. Once the algorithm converges, then each 𝜃
(5)
(for 𝑖 sufficiently large) will
be a sample from the posterior distribution. To that end, we discarded the first 9,000
iterations and only kept the remaining 1000 iterations to simulate the model. Given that this
approach to parameter estimation is probabilistic, we simulated the MH algorithm 200
independent times (for 10,000 iterations each) with a random initial guess.
2.3.4 Construction of magnitude of network activation metric
We defined network activation to allow us to compare the magnitude of signaling
across the three pathways (CD16, 2B4 and NKG2D). While the individual phospho-species
are known to contribute to specific cellular functions involved in NK cell signaling
3,34,46,54,58–
60
, the scope of our work is to compare the effectiveness of stimulating one pathway versus
another. Thus, we first determine which species are important to consider in terms of
activation of the network.
Based on literature evidence, we determined the following five species are crucial for
activating the NK cell based on experimental studies: (1) pErk, (2) pAkt, (3) pPLCγ, (4) pVav
and (5) pSLP76. The magnitude of network activation must relate to the magnitude of
activation of the above species. Hence, we concatenate the above species’ concentrations
over time into a vector:
g(𝑡) = xpErk(𝑡),pAkt(𝑡),pPLCγ(𝑡),pVav(𝑡),pSLP76(𝑡)y.
23
Since the individual phospho-species’ concentrations are continuous with respect to time 𝑡,
g(𝑡) is also continuous with respect to 𝑡 and thus measurable
61
. By construction, the arc g is
a function that maps the time interval (in minutes) [0,60] into ℝ
<
; that is to say, g : [0,60] →
ℝ
<
. We used the Bochner-norm to define the magnitude of g(𝑡). Briefly, the Bochner-norm
of is defined by:
‖g(𝑡)‖
=
+
(1,?1)
≝‖g(𝑡)‖
ℝ
,
A
?1
1
d𝑡
2
A
B
.
Since the image of g is an element of ℝ
<
, and since all norms define on ℝ
C
form
equivalence classes
62
, we set 𝑝 =1 and use the 𝐿
2
norm on ℝ
<
for simplicity. In addition, each
component 𝑔
5
is non-negative, has finite measure and the sum of the components of g, for all
𝑡 ∈ [0,60], is finite. Thus, the Bochner-norm simplifies to
‖g(𝑡)‖ = ‖g(𝑡)‖
ℝ
, d𝑡
?1
1
= ¡ 𝑔
5
(𝑡) d𝑡
?1
1
<
562
.
We arrive at the above metric for the magnitude of network activation, which is simply the
sum of the magnitude of activation of the individual phospho-species, as given by the area
under the curve for the species’ concentration profile. We used the MATLAB function trapz
(which uses trapezoidal numerical integration) to estimate the area under the curve for each
component of g(𝑡). We acknowledge that how much each phospho-species contributes to NK
cell activation is unknown. For simplicity, we assume equal contribution from the individual
24
species. This assumption can be adjusted as data for the contributions of each species
become available.
2.3.5 Clustering and principal component analysis
We used the built-in MATLAB functions kmeans and pca to perform k-means
clustering and principal component analysis, respectively. Briefly, k-means clustering
63
allows us to partition a given dataset into k clusters using the (default) Euclidean distance
metric. Principal component analysis
64
enables us to project a given dataset on to a new
coordinate system where each coordinate is a linear combination of the original variables in
the dataset. Moreover, the principal components (i.e., new coordinates) are selected such
that they maximize the total variance in the data. These approaches are used to determine
which estimated parameter sets are similar to one another.
2.4 Results
2.4.1 Model of NK cell signaling matches experimental data
We generated a mathematical model of NK cell signaling that includes three main
pathways: CD16, 2B4 and NKG2D. When these receptors are stimulated, they activate the
cell via cascades of phosphorylation reactions (Figure 2-1): activation of the SFK, facilitated
by the ligand-bound phosphorylated receptors, catalyzes the activation of the Akt, SLP76-
Vav-Erk and PLCγ pathways. We simulated these reactions in the form of nonlinear ODEs
25
using established Michaelis–Menten kinetics. The model is provided in Supplementary File
S1. The model was calibrated to immunoblot data
26,32–34
, where we quantified the temporal
change in the optical density of protein bands from images of immunoblot experiments using
ImageJ
55
. Specifically, we used the normalized levels of the following phosphorylated
species: pSFK, pZAP70, pLAT, ppSLP76, pPLCγ, pVav, pErk, pAkt and SLP76 phosphorylated
at Y113 and Y128. We calibrated the model predictions by estimating the parameter values
using a Bayesian perspective
56
, and by implementing the Metropolis–Hastings algorithm
(see section 2.3.3). In brief, the model parameters (83 in total) were estimated 200 times
using randomized initial guesses by fitting to experimental data. Moreover, the model
predictions were validated using a separate dataset. The combined error for each run can be
found in Supplementary Figure S1 (Figure S2-1). We proceeded with the 14 best parameter
sets that provided the lowest total error and simulated the model using these sets.
Interestingly, in initial simulations, we found that the 14 parameter sets led to
different responses with respect to network activation. To determine the dominant behavior
generated by the model, we clustered the network activation predicted by the 14 sets using
the kmeans and pca functions in MATLAB (see section 2.3.5). Our results are shown in
Supplementary Figure S2 (Figure S2-2), where we identified three unique clusters that
correspond to the degree of species activation (i.e., response) predicted by the model. To
ensure that the model predictions agree with experimental observations, we discarded the
parameter sets that yielded predictions inconsistent with NK cell signaling and cytotoxicity
studies
26,32–35,54,59,65
. Specifically, we removed parameter sets that induced a low amount of
species activation (i.e., < 1% of the species’ initial concentration was activated) and that did
not show a dose-dependent response when the ligand concentrations were changed.
26
This refined the 14 parameter sets down to five, which are found in the medium
response cluster in Supplementary Figure S2 (Figure S2-2). Specifically, the parameter sets
in this cluster are parameter sets 3, 5, 6, 12 and 13 in Supplementary Figure S1 (Figure S2-
1). The parameter distributions for the best set (i.e., lowest total error) of the five (parameter
set 3 in Supplementary Figure S1 (Figure S2-1)) are shown in Supplementary Figure S3
(Figure S2-3) using the final 1,000 iterations, illustrating that the parameters are well
behaved: the distributions are unimodal, and the values lie within a tight range. We used the
last 1,000 iterations from parameter estimation to simulate the model. In addition, the
posterior distribution of the parameters can be found plotted with the prior distribution in
Supplementary Figure S4 (Figure S2-4) – Supplementary Figure S7 (Figure S2-7). Since the
Metropolis–Hastings algorithm is a Monte Carlo Markov Chain (MCMC) algorithm, we
provide diagnostic information in the form of trace plots to demonstrate that the parameters
are identifiable. Supplementary Figures S8 (Figure S2-8) and Supplementary Figure S9
(Figure S2-9) show the traces for the best parameter set, and the traces for the remaining
four sets in the cluster are shown in Supplementary Figure S10 (Figure S2-10) and
Supplementary Figure S11 (Figure S2-11). Collectively, the parameters in each set (i.e.,
independent MCMCs) do converge, but at different iterations. There are parameters that
converge as early as the 2000th iteration, while most converge closer to the 5,000
th
iteration,
and only very few require > 5,000 iterations to converge. The dashed lines in Supplementary
Figures S8 (Figure S2-8) and Supplementary Figure S9 (Figure S2-9) indicate the cutoff
(9,000
th
iteration). The left side is considered the burning-in phase, and those values were
discarded, while the values on the right side were kept for model simulation.
27
The simulated concentration profiles are consistent with the training data (Figure 2-
2A – H). These results demonstrate the model predictions are in accord with the
experimental data for mono-stimulation of CD16 (blue lines), 2B4 (purple lines) and NKG2D
(orange lines). This is expected, since those data were used in model training to determine
the parameter values. To validate the model, we compared the model predictions to separate
experimental data not used during training. In particular, we quantified the optical density
of intracellular species from immunoblot images when 2B4 and NKG2D were simultaneously
stimulated at equal ligand concentrations. The results from model validation are shown in
Figure 2-2 | Model calibration and validation. The model was fit to experimental data for (A)
pSFK, (B) pZAP70, (C) pSLP76, (D) pLAT, (E) pPLCγ , (F) pVav, (G) pErk and (H) pAkt. The
modelpredictionswerevalidatedagainstseparatedatafor(I)pErk,(J)pAkt,(K)pPLCγand(L)
pVav under co-stimulation of 2B4 and NKG2D. Blue: CD16 pathway; purple: 2B4 pathway;
orange: NKG2D pathway; Green and Brown: 2B4 and NKG2D co-stimulation from separate
experiments. Note that the green and brown lines represent independent western blot
experiments that only differ in the time-points of data collection. Solid lines: mean model
predictions from 1000 parameter estimates. Shaded area: standard deviation of mean model
predictions.Squares,circlesandtriangles:experimentaldata.
28
Figure 2-2I – L. The model captures the signaling dynamics of several species upon co-
stimulation of 2B4 and NKG2D.
We note that to perform parameter estimation, we separated the time-course
datasets into two groups based on which datasets involved mono-stimulation of the
receptors and which involved co-stimulation. Currently, we have 21 datasets: 15 involving
mono-stimulation and the remaining six sets involving co-stimulation. We designated the
mono-stimulation datasets to be used for model calibration, keeping the co-stimulation
datasets for model validation. Although the results from parameter estimation can be
sensitive to which datasets are used for model training versus model testing, the total
number of combinations of selecting 15 datasets for training from full set of 21 datasets is
quite large r
72!
20 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 𝛼 and 𝛽 affect the shape and scale of the noise,
respectively. Also, note that the prior distributions do not affect the acceptance ratio since
APQ
(!)
R
A(Q
∗
)
=1 regardless of the parameter vectors 𝜃
(O)
or 𝜃
∗
and independent of the arbitrary
intervals [𝑎
5
,𝑏
5
] for each 𝑖
$%
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 𝜃
(O)
and 𝜃
∗
. 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 2,000
th
and 5,000
th
iteration. Given that the MH algorithm is a
stochastic optimization method, we simulated the MH algorithm 200 times using
independent random initial guesses. In addition, the uncertainty in the parameter estimates
from our previous model
86
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
1,000 iterations of the MH algorithm, which we take to be samples from the posterior
distribution. The model is provided in supplementary material file S1, and the list of model
species, reactions, and parameters is provided in supplementary material file S2.
72
3.3.4 Information-theoretic sensitivity analysis
We employed an entropy-based sensitivity analysis for informing which parameters
(model inputs) share a significant degree of mutual information with the amount of secreted
GZMB and PRF1 (model outputs). We follow the methods described previously by Lüdtke et
al
98
. Entropy, in the sense of Shannon, is described as the average information content of a
random variable. That is, for any random variable 𝑌, the (Shannon) entropy of 𝑌 is given by
𝐻(𝑌)=−¡𝑝(𝑦)log
7
x𝑝(𝑦)y
S
.
The conditional entropy of 𝑌 given 𝑋 =𝑥 x𝐻(𝑌|𝑋 =𝑥)y is defined analogously by
using the conditional probability 𝑝(𝑦|𝑋 =𝑥). Moreover, the quantity 𝐻(𝑌|𝑋) measures the
average amount of information remaining 𝑌 given that we observed another random
variable 𝑋. Specifically,
𝐻(𝑌|𝑋)=¡𝑝(𝑥)𝐻(𝑌|𝑋 =𝑥)
T
=−¡𝑝(𝑥,𝑦)log
7
°
𝑝(𝑥,𝑦)
𝑝(𝑥)
±
T,S
.
Note that 𝐻(𝑌)≥𝐻(𝑌|𝑋) as any random variable 𝑋 can only explain away some of
the information in 𝑌, if any at all. If 𝐻(𝑌|𝑋)=𝐻(𝑌), then 𝑌 is independent of 𝑋. Alternatively,
if 𝐻(𝑌|𝑋) =0, then knowing 𝑋 completely determines 𝑌. As defined in the study by Lüdtke,
et al.
98
, conditional entropies of the form 𝐻(𝑌|{𝑋
2
,…,𝑋
C
} \ 𝑋
5
) measure the total effect that
73
a particular input 𝑋
5
exerts on the output 𝑌. Furthermore, the quantity 𝐻(𝑌|{𝑋
2
,…,𝑋
C
})
determines the amount of information remaining in 𝑌 once we observed all the inputs. This
can be thought of as the residual information that persists in 𝑌 for which the inputs 𝑋
2
,…,𝑋
C
cannot account for. Thus, the total order sensitivity index for each input 𝑋
5
is defined by
𝑆
5
=
𝐻(𝑌|{𝑋
2
,…,𝑋
C
} \ 𝑋
5
)
𝐻(𝑌)−𝐻(𝑌|{𝑋
2
,…,𝑋
C
})
.
Indeed, inputs with a higher sensitivity index suggest that the output is sensitive to
variations in the input. In our case, 𝑋
5
represents the kinetic parameters in our model,
whereas 𝑌
2
and 𝑌
7
represent the amount of secreted GZMB and PRF1, respectively, after 60
minutes of receptor stimulation to mimic the experimental conditions from the study by
Srpan et al
53
. We drew 250 random samples for each 𝑋
5
(independently) using a uniform
distribution on the interval [0.5E(𝑋
5
),1.5E(𝑋
5
)], where E is the expectation operator and the
distribution of each 𝑋
5
is given by parameter estimation (see section 3.3.3) or from the study
by Makaryan and Finley
86
. Next, we simulated the model using these 250 samples to generate
a distribution for each of the outputs 𝑌
2
and 𝑌
7
to then compute the total order sensitivity
indices.
3.3.5 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.
29
engineered novel genetic circuits in
74
mammalian cells to respond to inflammatory signals (i.e., IL-1𝛽) by producing anti-
inflammatory proteins. In addition, the use of synthetic biology methods with clinical
applications is the subject of recent reviews
99–101
. The synthetic Notch (synNotch) signaling
pathway
102
, 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
71,75,83
. 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. 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.
75
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
103–106
. Recently, lncRNA
pulldown assays
106–108
, 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 protein sequestering lncRNAs is that, once discovered, they can be sequenced
and reverse-transcribed via polymerase chain reaction (PCR) to manufacture their
complementary DNA strand (cDNA)
109
. 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 (𝜇M×min)
J2
and (2)
there are two plasmids (at fixed amounts) controlling the expression of lncRNA and the
76
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
𝐾
-
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
110
and RNA degradation
111,112
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.
3.3.6 Optimization of 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,
77
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
62,113–115
; specifically, we minimized the following stochastic objective
function, given the model parameters 𝜃:
min
U,V,W
0
·(1−𝜌)(u+v+𝑅
1
)−𝜌º¡ExGZMB
X
(60,𝜃)y+ExPRF1
X
(60,𝜃)y
Y
X62
»¼.
Indeed, we minimize the objective function using the sample average approach
(SAA)
116–119
. Here, u, v, and 𝑅
1
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≤1,000 (copies×cell
J2
)
0≤𝑅
1
≤10 𝜇M
. We set the
upper bound of the plasmid concentrations based on what is defined as a high copy
number
29,101,102
. In addition, the upper bound on the initial value of synNotch receptor is
based on a study
120
where Chinese hamster ovary (CHO) cells were genetically modified to
maximally produce human IgG; the values ranged from 0.3 – 20 𝜇M, from which we chose 10
𝜇M. 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
78
of 𝑅
1
was either 20 𝜇M 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 minutes of receptor stimulation, where 𝑟 =1,…,𝑁 is the number of
rounds of receptor stimulation and the expectation is taken over the parameters 𝜃. Since
GZMB and PRF1 are non-negative, the second term is, in fact, a maximization problem given
that argmin
−𝑓 =argmax
𝑓 for all non-negative 𝑓. The constant 𝜌 ∈(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 ri.e., 𝜌 =
2
7
t.
To solve this optimization problem, we used the mesh-adaptive directed 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.
79
3.4 Results
3.4.1 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
86
. 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 (Figure 3-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 PLC𝛾 pathways. In particular, we
assume that pVav and pPLC𝛾 mediate the secretion of GZMB and PRF1, as they have been
correlated with NK cell cytotoxicity
39
and a discharge of intracellular calcium ions
1,26,28,35,95
,
which mediate exocytosis, respectively. Although Akt and Erk are necessary for cell survival
and proliferation
26,28,34,35
, respectively, their contribution to NK cell secretion of GZMB and
PRF1 remains elusive. As it stands, given the available data, we consider pVav and pPLC𝛾 as
the mediators for NK cell secretion.
The model was calibrated to data from the study by Srpan et al.
53
using a Bayesian
perspective to parameter estimation
56,57
; namely, we implemented the Metropolis–Hastings
80
(MH) algorithm (see section 3.3.3) 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 Figure S3-1.
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 Figure S3-1). The trace plots for each parameter in the best run are shown in Figure S3-2,
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 1,000 times
using the final 1,000 iterations of each parameter from the best run (Figure S3-3) 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 (Figure 3-2). The model can reproduce the time evolution of secreted PRF1
mediated by stimulation of NKG2D (Figure 3-2A) and CD16 (Figure 3-2B). In addition, the
model can replicate sequential stimulation data where NK cells are stimulated via one
pathway for two consecutive rounds, each for 60 minutes, and then stimulated via the other
pathway for the third round of stimulation. Under these conditions, the simulated
concentrations of secreted PRF1 (Figure 3-2C and Figure 3-2D), the receptors (Figure 3-2E
– H), and intracellular PRF1 (Figure 3-2I – J) all agree with the experimental data. To further
validate the model predictions, we simulated the concentration of secreted PRF1 when both
Rituximab and MICA were decreased to 1 𝜇g/ml (Figure 3-2K), and these model simulations
81
agree with experimental measurements as well. In addition, we demonstrate in Figure S3-4
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
Figure 3-2 | Model training and validation. The model was trained to, and tested against, in vitro
NK cell stimulation data from the study by Srpan et al. In all panels, blue markers
and bars represent mean experimental data, while red and green signify the mean model
predictions for the NKG2D and CD16 pathways, respectively. Error bars indicate one standard
deviation. Except for panel (k), the ligand concentrations used to stimulate NKG2D and CD16 are
2.5 𝜇g/ml of MICA and 10 𝜇g/ml of Rituximab, respectively. (a) and (b) Normalized time series
data of PRF1 secretion using (a) 2.5 𝜇g/ml of MICA or (b) 10 𝜇g/ml of Rituximab. (c) and (d)
Normalized PRF1 secretion per round of stimulation (60 minutes each) by stimulating the CD16
or NKG2D pathway and then stimulating the NKG2D pathway or CD16 pathway. (e) and (f)
Normalized concentration of CD16 per round of stimulation (60 min each) by stimulating the
CD16 or NKG2D pathway and then stimulating the NKG2D or CD16 pathway. (g) and (h)
Normalized concentration of NKG2D per round of stimulation (60 min each) by stimulating the
CD16 or NKG2D pathway and then stimulating the NKG2D or CD16 pathway. (i) and (j)
Normalized concentration of intracellular PRF1 per round of stimulation (60 min each) by
stimulating the NKG2D or CD16 pathway and then stimulating the CD16 or NKG2D pathway. (k)
Normalized concentration of PRF1 after 60 min of stimulation using different concentrations of
ligand.
82
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 (Figure S3-5) indicate that the best fit (Run 1 from Figure S3-1) 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 Figure S3-5 generated prediction
errors similar in magnitude to the original combination. We found that the parameter
distributions (from the final 1,000 iterations of the MH algorithm) from combinations 4, 8,
and 12 were comparable to the original combination with a significant overlap (Figure S3-
6). 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 Figure
S3-3, which overlap with the distributions in Figure S3-6. Overall, we generated an
experimentally validated mechanistic model of NK cell degranulation, which can recreate
results from the study by Srpan et al.
53
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.
3.4.2 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 section 3.3.4). This technique, as shown by Lüdtke et
83
al.
98
, 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. 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 (Figure S3-7). The fifteen most influential
parameters are found in the pSFK-pSHP-pVav-pPLC𝛾 subgraph, where each parameter can
explain more than 35% of the information in GZMB secretion promoted by the CD16
receptor (Figure 3-3A; Figure S3-7A). This also holds true for NKG2D-mediated GZMB
secretion (Figure 3-3B; Figure S3-7C) as well as for PRF1 secretion promoted by either
receptor (Figure S3-7B and Figure S3-7D). These results suggest that this subnetwork
(Figure 3-3C) 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
Figure 3-3C. Thus, we intervened on these parameters (individually) by modifying their
84
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 pPLC𝛾), 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 minutes of receptor stimulation (Figure 3-3D and Figure S3-8). The simulated percent
Figure 3-3 | Sensitivity analysis shows that inhibition of pSHP activation is most influential for
GZMB and PRF1 secretion. Total order sensitivity indices of the top fifteen influential parameters
for the amount of GZMB secretion after 60 min of stimulation of the (a) CD16 pathway or (b)
NKG2D pathway. (c) Almost all the influential parameters from the sensitivity analysis can be
found in the pSFK-pSHP-pPLC𝛾-pVav sub-network of the full model. (d) Predicted percent change
in GZMB secretion when decreasing the catalytic rate constants 𝑘
123
, parameters involved in the
sub-network depicted in (c) after 60min of stimulation of CD16 (left) and NKG2D (right). In panel
(d): green, decreasing pSHP activation; blue, decreasing pSFK deactivation; magenta, decreasing
pVav deactivation; cyan, decreasing pPLC𝛾 deactivation; markers, mean model prediction; and
error bars, one standard deviation.
85
increase in GZMB secretion is greatest when the pSFK → pSHP reaction is inhibited (Figure
3-3C; decreasing 1) and when either CD16 or NKG2D (Figure 3-3D) is stimulated. The model
predicts that reducing the catalytic rate constant for this pathway has a larger effect when
CD16 is stimulated, compared to NKG2D stimulation (Figure 3-3D). The remaining
interventions were not as influential, however (Figure S3-8A and Figure S3-8D).
Furthermore, this conclusion is equally true for PRF1 secretion (Figure S3-8B, Figure S3-8C,
Figure S3-8E, and Figure S3-8F). These results suggest that disrupting the incoherent feed-
forward loop (IFFL) pSFK → pSHP ⊣ pX, where X is either Vav or PLC𝛾, 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 Figure 3-3C.
3.4.3 The optimal synNotch system depends on the number of rounds of
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 (Figure S3-9).
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
41,110
. Ultimately, this causes NK cells to become less cytotoxic over time.
Therefore, to maximize GZMB and PRF1 secretion over multiple rounds of stimulation, we
86
need to induce the production of these molecules, in addition to reducing SHP activation.
Fortunately, the field of synthetic biology
99–101
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
102
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
99–101,121
that can be tailored to a specific stimulus
(Figure 3-4A). The NK cell signaling network promoted by the endogenous receptors CD16
and NKG2D combined with the synthetic Notch system is illustrated in Figure 3-4B. 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
103,104,106
(see section
3.3.5). 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.
87
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
Figure 3-4 | Optimal strategies depend on the number of rounds of stimulation. (a) Diagram of the
synNotch signaling pathway. 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.
88
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
29,101,120
, 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 section 3.3.6).
Interestingly, the optimized synNotch system has different characteristics when
paired with CD16 (Figure 3-4C) or NKG2D (Figure 3-4F). When CD16 is stimulated for 1–2
rounds only (Figure 3-4C), 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 hours) 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 (Figure 3-4F). In contrast, when CD16 is
stimulated for 3 – 5 rounds, the optimal amounts for each plasmid and 𝑅
1
(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
89
molecule-coding plasmid and 𝑅
1
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 velocity of phospho-
protein decay
26,122,123
(see Table I, 𝑘
FGH
). 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 (Figure 3-4F) –
similar to CD16 (Figure 3-4C). 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 𝑅
1
does not immediately reach its maximum value. This difference between
the optimal amounts of 𝑅
1
is due to the initial amounts of CD16 and NKG2D receptors (38
and 0.3 𝜇M, respectively), compared to the maximum value that 𝑅
1
can have, 10 𝜇M. 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 𝑅
1
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 𝑅
1
, a large 𝑅
1
would lead to less secretion by directing the
90
input signal more toward cytolytic molecule production. Moreover, since the optimal 𝑅
1
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 𝑅
1
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.
3.4.4 The optimal CD16- and NKG2D-synNotch pairs show qualitative
differences
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
91
“semi-max control,” respectively. For CD16, we used the optimal amounts from Figure 3-4C
for rounds 1, 3, and 6, respectively, for these three conditions, whereas for NKG2D, we used
the optimal amounts from Figure 3-4F for rounds 1, 11, and 15, respectively.
When considering GZMB secretion (Figure 3-4D), 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
(Figure 3-4G). This counterintuitive result that no control is optimal is due to the
competition for ligand between the endogenous and the synNotch receptors as well as the
time delay in protein production. For short-term stimulation, it is best for the NK cell to rely
on endogenous signaling to maximize GZMB and PRF1 secretion (Figure 3-4D and Figure 3-
4G). These simulated results help contextualize why it was optimal to not intervene in the
first few rounds of stimulation (Figure 3-4C and Figure 3-4F).
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
53
. 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 𝑘
FGH),( [\]7E
(Figure S3-3), which is approximately three times
92
larger than 𝑘
FGH),( &E2?
, 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 (rounds 3 and 11 for CD16 and NKG2D
stimulation from Figure 3-4C and Figure 3-4F, respectively) leads to more GZMB secretion
for 2 – 4 rounds of CD16 stimulation and for 3 – 10 rounds of NKG2D stimulation (green
squares in Figure 3-4D and Figure 3-4G) when compared to the no control set (gray squares
in Figure 3-4D and Figure 3-4G). Finally, the semi-max control set (rounds 6 and 15 for
stimulation of CD16 and NKG2D from Figure 3-4C and Figure 3-4G, respectively) secretes
more GZMB after 5 rounds of CD16 stimulation and after 11 rounds of NKG2D stimulation
when compared to the max control set (compare green with orange squares in Figure 3-4D
and Figure 3-4G). As a result, the semi-max control set is predicted to be the preferred
strategy when the number of rounds of stimulation continues to increase. These conclusions
hold for PRF1 secretion also (Figure 3-4E and Figure 3-4H).
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 (Figure
3-4D). 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
93
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 𝑅
1
is greater than the initial value of NKG2D (0.3 𝜇M), 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
𝜇M), in contrast, is greater than 𝑅
1
, 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 (Figure 3-4D, Figure 3-
4E, Figure 3-4G, and Figure 3-4H), the secreted amount of cytolytic molecules decreases. We
interrogated the model to better understand this prediction. While the intracellular pool of
GZMB (Figure S3-10A) and PRF1 (Figure S3-10B) continues to increase as we increase the
number of rounds of stimulation, the intracellular amount of pPLC𝛾 (Figure S3-10C) and
pVav (Figure S3-10D) eventually decreases to zero. Since the secretion of GZMB and PRF1 is
mediated by pPLC𝛾 and pVav in our model, it follows that once the amount of pPLC𝛾 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.
94
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 (Figure S3-11) or NKG2D (Figure S3-12). 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 Figure 3-4C and Figure 3-4F. 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 (Figure S3-11A). Moreover, when
compared to the parameters in the endogenous pathway (Figure S3-11A; 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 Figure 3-3A, which belong to
the endogenous pathway. Similarly, these conclusions hold for the semi-max control case of
GZMB secretion (Figure S3-11C) as well as for PRF1 secretion (Figure S3-11B and Figure S3-
11D) and GZMB and PRF1 secretion promoted by the NKG2D pathway (Figure S3-12). To
further verify these results, we intervened on the 𝑅
1
and 𝑘
'(
parameters in the synNotch
pathway (which are the most influential parameters from the synNotch system) as we
presented above (Figure 3-3D) and compared the results with the decreasing activation rate
of pSHP (Figure S3-13). We found no significant change in the secretion of GZMB or PRF1. In
95
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.
3.4.5 Transcription factor cooperativity impacts the performance of synNotch
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 Figure 3-4C and Figure 3-4F 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
124,125
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):
rate of RNA production ≝𝑟(𝑡)=𝑘
$)"$
u
2
2/L
1
2
34(5)
N
*
, where u is the amount of plasmid and 𝑘
$)"$
is the transcription rate. By taking u and 𝑘
$)"$
to be fixed, we observe that the rate of RNA
production increases as we decrease the parameters n (the Hill coefficient) and 𝐾
-
(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
126–128
(i.e., the TF). It is possible to change the rate of RNA production by varying n. For
96
example, decreasing n should increase the velocity of RNA production. Indeed, if we take, for
example, n
1
≝0.5n and n
2
≝1.5n, then,
𝑟
1
(𝑡) =𝑘
$)"$
u
1
1+¤
𝐾
-
TF(𝑡)
¥
(
0
=𝑘
$)"$
u
1
1+¤¤
𝐾
-
TF(𝑡)
¥
(
¥
1.<
>𝑟(𝑡),
where 𝑟(𝑡) is the rate of RNA production when n = 1, and similarly,
𝑟
2
(𝑡)=𝑘
$)"$
u
1
1+¤
𝐾
-
TF(𝑡)
¥
(
6
<𝑟(𝑡) <𝑟
1
(𝑡),
for all time 𝑡. 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 (Figure 3-5A) and cytolytic molecule-coding plasmid (Figure 3-5B) are at their
maximal values and 𝑅
1
(Figure 3-5C) is near-maximal even with one round (i.e., 1 hour) 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
97
(n<1) with darker shade (n>1) in Figure 3-5A – C). These results hold for the NKG2D
pathway as well (Figure S3-14A – 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 Figure 3-5A – 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.
Figure 3-5 | Transcription factor cooperativity impacts CD16-synNotch performance. (a)–(c)
Optimal values of (a) lncSHP-coding plasmid, (b) cytolytic molecule-coding plasmid, and (c) initial
synNotch for negative cooperativity (lightest shade), no cooperativity (intermediate shade), and
positive cooperativity (darkest shade). (d)–(f) The amount of secreted GZMB from CD16-
synNotch pair when using the optimal values from (d) round 1, (e) round 3, and (f) round 6 from
(a)–(c). (g)–(i) The amount of secreted PRF1 from CD16-synNotch pair when using the optimal
values from (g) round 1, (h) round 3, and (i) round 6 from (a)–(c). (d)–(i) Markers, mean model
prediction; error bars, one standard deviation; lightest shade, negative cooperativity (n=0.5);
intermediate shade, no cooperativity (n=1); darkest shade, positive cooperativity (n=1.5).
98
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 Figure 3-4C. We apply the model by using the
predicted optimal values from the same rounds with the different Hill coefficients. When we
consider GZMB secretion (Figure 3-5D – 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 (Figure 3-5G – I) and the NKG2D pathway (Figure S3-
14D – 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 Figure S3-14E and Figure S3-14F with Figure 3-5E and Figure 3-5F),
suggesting that increasing RNA production shifts the response from small-but-sustained to
large-and-transient (compare the lightest shade with intermediate shade in Figure S3-14E,
Figure S3-14F, Figure S3-14H, and Figure S3-14I). In conclusion, we found that TF negative
cooperativity is predicted to augment the performance of synNotch, whereas positive
cooperativity is predicted to reduce it.
99
3.4.6 The synthetic-endogenous pathway pairs promote a biphasic 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 𝜇M and simulate the model for
one round of stimulation until the steady state is attained. That is, instead of stimulating the
receptors for 𝑘 rounds (60 minutes each) with 𝑘−1 intermediary washing steps (15
minutes 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 𝑘
FGH
, 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
(Figure 3-6A) 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 𝜇M—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
100
of GZMB (Figure 3-6A), the benefits of synNotch are observed when Rtx is between 10
-3
and
10 𝜇M and peak at approximately 0.7 and 0.2 𝜇M 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 (Figure 3-6B) or via the NKG2D-synNotch
pair (Figure 3-6C – D).
Next, we probed the model to better understand the cause of the biphasic, steady state
response. In Figure S3-15A – D, we demonstrate the time evolution of secreted GZMB via the
Figure 3-6 | Long-term behavior of the combined synthetic-endogenous pathway. (a) and (b) The
steady state levels of (a) GZMB and (b) PRF1 mediated by the CD16-synnotch combination as a
function of the stimulus. (c) and (d) The steady state levels of (c) GZMB and (d) PRF1 mediated
by the NKG2D-synNotch tandem as a function of the stimulus. Solid line, mean model prediction;
shaded area, one standard deviation; gray, no control set [round 1 for CD16 and NKG2D from
Figure 3-4C and Figure 3-4F, respectively]; green, max control set [round 3 for CD16 and round
11 for NKG2D from Figure 3-4C and Figure 3-4F, respectively]; orange, semi-max control set
[round 6 for CD16 and round 15 for NKG2D from Figs. 4(c) and 4(f), respectively].
101
CD16-synNotch pair at four distinct levels of Rtx: 10
-2
, 10
-2
, 1, and 10 𝜇M, respectively. Also,
we show the time evolution of secreted PRF1 in Figure S3-15E – H. We use the same control
sets as we defined above. In these figures, we observe that when Rtx is 10
-1
𝜇M, 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, Figure S3-15A with Figure S3-15B – 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 hours when
Rtx is equal to 10
-2
𝜇M (Figure S3-15I) and close to 8 hours when Rtx is equal to 10
-1
𝜇M
(Figure S3-15J). However, once the concentration of Rtx is 1 or 10 𝜇M (Figure S3-15K – I),
the concentration of phospho-CD16 undergoes a large rapid increase but decreases
precipitously after a 3- or 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 𝑘
_($ &E2?
, Table I) is proportional to the amount of
phospho-CD16. That is, degradation of phospho-CD16 ∝ 𝑘
_($ &E2?
pCD16(𝑡), where 𝑘
_($ &E2?
is between 1.4 and 1.5 min
-1
(Figure S3-3B). 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 (Figure S3-16). 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.
102
3.5 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
53
.
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
103
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 ⊣ pX, where X is either Vav or PLC𝛾. Interestingly, this motif is not
uncommon in biological networks
129–132
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
69,70,133–135
. 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 plasmid encoding lncRNA for the phosphatase, the
104
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. 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
136–138
between the TF and its DNA
binding site on the plasmid (i.e., decreasing 𝐾
-
).
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
𝜇M or
greater than 10 𝜇M, the added benefits of synNotch are nullified. In fact, the speed of
degradation of the intracellular signaling molecules is proportional to the amount of
105
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
77
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
53,77,139–141
, 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
106
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. 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.
53
, 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 𝜌 determined how much emphasis we placed
107
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
108
provides actionable insight into engineering robust NK cells with applications to
immunotherapies.
3.6 Acknowledgements
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).
3.7 Appendix
109
Figure S3-1 | Error in model training and testing. Model parameters were estimated via the
Metropolis-Hastings (MH) algorithm for 200 independent runs (for 10,000 iterations each),
shown as the horizontal axis. The final 1,000 iterations were used to simulate the model and
compare the results with the experimental data. We took the squared Euclidean norm to
estimatethediscrepancyanddividedbythenumberofdatapoints(N).Notethatthenumberof
datapointsisdifferentformodeltraining(N=34)andtesting(N=16).Bars,meanvalue;error
bars, one standard deviation. The 200 independent runs are sorted based on the (mean) total
error.
110
Figure S3-2 | Trace plots of estimated parameters. A scatter plot of the estimated parameters
depictingtheparameter’svalueateachiteration.Theseresultsaretakenfromthebestfit(Run1
fromFigureS1)fromthe200independentsimulationsoftheMetropolis-Hastingsalgorithm.
111
FigureS3-3|Histogramsofestimatedparameters.Theestimatedposteriorprobability
distributionoftheparametersconditionalontheexperimentaldata.Theseresultsshowthefinal
1,000iterationsoftheMetropolis-Hastingsalgorithmfromthebestfit(Run1fromFigureS1).
112
FigureS3-4|Modeltrajectoriesincomparisonwithexperimentalconditions.SimilartoFigure2,
we show the time series projections of the model against the experimental data. In all panels,
bluemarkersandbarsrepresentmean experimentaldata,while redandgreen signify themean
model predictions for the NKG2D and CD16 pathways, respectively. Error bars indicate one
standard deviation. Shaded area of model predictions represents one standard deviation. Gray
areas show intermediary washing step (15 minutes each). Except for panel (K), the ligand
concentrationofMICAandRituximab(Rtx)are2.5Qg/mLand10Qg/mL,respectively.(AandB)
Normalized time series data of PRF1 secretion using (A) 2.5 Qg/mL of MICA or (B) 10 Qg/mL of
Rituximab.(CandD)NormalizedPRF1secretionperroundofstimulation(60minuteseach)by
stimulating the CD16 or NKG2D pathway and then stimulating the NKG2D pathway or CD16
pathway. (E and F) Normalized concentration of CD16 per round of stimulation (60 minutes
each) by stimulating the CD16 or NKG2D pathway and then stimulating the NKG2D or CD16
pathway. (G and H) Normalized concentration of NKG2D per round of stimulation (60 minutes
each) by stimulating the CD16 or NKG2D pathway and then stimulating the NKG2D or CD16
pathway. (I and J) Normalized concentration of intracellular PRF1 per round of stimulation (60
minutes each) by stimulating the NKG2D or CD16 pathway and then stimulating the CD16 or
NKG2D pathway. (K) Normalized concentration of PRF1 after 60 minutes of stimulation using
differentconcentrationsofligand.
113
Figure S3-5 | Randomized model calibration. The data used for model training and validation
was randomized into sixteen different combinations including the original case (Figure S1). For
each combination, the model parameters were estimated via the Metropolis-Hastings (MH)
algorithmfor8independentruns(for10,000iterationseach),shownasthehorizontalaxis.The
final 1,000 iterations were used to simulate the model and compare the results with the
experimentaldata.WetookthesquaredEuclideannormtoestimatethediscrepancyanddivided
by the number of data points (N). Note that the number of data points is different for model
training(N=34)andtesting(N=16).Bars,meanvalue;errorbars,onestandarddeviation.
114
Figure S3-6 | Histograms of parameter estimates from combinations with lowest total error
overlap considerably. The estimated posterior probability distribution of the parameters
conditional on the experimental data. These results show the final 1,000 iterations of the
Metropolis-Hastings algorithm from the original combination, and combinations 4, 8 and 12
fromFigureS5.Thesedatashowthattheparametersarestablewithrespecttovariationsinthe
modelcalibrationprocess.Theparametervaluesusedtosimulatethemodel(FigureS3)overlap
withtheestimatesaswell.
115
FigureS3-7|Totalordersensitivityindicesofkineticparameters.AsinFigure3A-B,we
show the total sensitivity index of all 89 parameters in the model with respect to the amount of
(AandC)GZMBand(BandD)PRF1secretionafter60minutesofstimulationof(AandB)CD16
and(CandD)NKG2D.
116
Figure S3-8 | Intervention on remaining parameters found in Figure 3A, B. Based on the results
from Figure 3, we found 8 different reactions that can lead to an increase the amount of GZMB
and PRF1 degranulation after 60 minutes of receptor stimulation. We varied the catalytic rate
constants for the significantly influential parameters identified by sensitivity analysis and
quantified the percent change in (A) GZMB and (B – C) PRF1 after 60 minutes of stimulation of
CD16and(D)GZMBand(E–F)NKG2Dafter 60minutesof stimulation of NKG2D. In all panels,
markersrepresentthemeanmodelpredictionanderrorbarsrepresentonestandarddeviation.
117
FigureS3-9|InhibitionofpSHPactivationdepletesintracellularpoolofcytolyticmolecules.The
simulated change in the absolute, intracellular concentration of GZMB (A) and PRF1 (B) after
receptor stimulation for 60 minutes. The horizontal axis depicts the percent decrease in the
catalytic rate parameter regulating SHP phosphorylation from baseline. Markers represent the
mean model prediction and error bars represent one standard deviation. Green, CD16
stimulation;Red,NKG2Dstimulation.
118
FigureS3-10|DepletionofintracellularpPLC9andpVavleadstoadecreaseofcytolytic
molecule secretion. The performance of the NK cell with and without controllers. The optimal
amount of plasmid and synNotch receptor are as in Figure 4. The simulated intracellular
concentration of (A) GZMB, (B) PRF1, (C) pPLCP and (D) pVav via the CD16 and NKG2D
pathways.Markersrepresentthemeanmodelpredictionanderrorbarsrepresentonestandard
deviation.Gray,“nocontrol”;Green,“maxcontrol”;Yellow,“semi-maxcontrol”(seemaintextfor
definitions).
119
FigureS3-11|SensitivityanalysisofsynNotchspecificparametersviaCD16pathway.The
totalordersensitivityindiceswhenusingthe“maxcontrol”valuesoptimizedforanintermediate
number of rounds of stimulation from Figure 4D (GZMB secretion) and Figure 4E (PRF1
secretion) for synNotch pathway parameters (left) and all model parameters (right) for (A)
GZMB secretion and (B) PRF1 secretion. Total order sensitivity indices when using the “semi-
max control” values optimized for long-term stimulation for synNotch pathway parameters
(left) and all model parameters (right) for (C) GZMB secretion and (D) PRF1 secretion. Orange
barsarespecifictothesynNotchpathway.
120
Figure S3-12 | Sensitivity analysis of synNotch specific parameters via NKG2D pathway. The
total order sensitivity indices when using the “max control” values optimized for an
intermediate number of rounds of stimulation from Figure 4D (GZMB secretion) and Figure 4E
(PRF1 secretion) for synNotch pathway parameters (left) and all model parameters (right) for
(A) GZMB secretion and (B) PRF1 secretion. Total order sensitivity indices when using the
“semi-max control” values optimized for long-term stimulation for synNotch pathway
parameters (left) and all model parameters (right) for (C) GZMB secretion and (D) PRF1
secretion.OrangebarsarespecifictothesynNotchpathway.
121
FigureS3-13|ComparisonofincreasingsynNotchreceptor( :
!
)anditsaffinitytoligand
( >
"#
)todecreasingpSHPactivation.TheoptimalinitialsynNotchreceptorconcentration( :
!
)
found in Figures 4B, C, along with the >
"#
parameter, were varied as in Figure 3. The percent
change in (A) GZMB secretion and (B) PRF1 secretion from baseline using the “max control”.
Thiscorrespondsto3and11roundsforCD16andNKG2Dstimulation,respectively.Thepercent
change in (C) GZMB secretion and (D) PRF1 secretion from baseline using the “semi-max
control”. This corresponds to 6 and 15 rounds for CD16 and NKG2D stimulation, respectively
Markers represent the mean model prediction and error bars represent one standard deviation.
Black,inhibitingpSHPactivation;red,increasing :
!
;blue,increasing >
"#
.
122
FigureS3-14|TranscriptionfactorcooperativityimpactsNKG2D-synNotchperformance.(A–C)
Optimalvaluesof(A)lncSHP-codingplasmid,(B)cytolyticmolecule-codingplasmid,and(C)and
initial synNotch with negative cooperativity (lightest shade), no cooperativity (intermediate
shade), and positive cooperativity (darkest shade). (D – F) The amount of secreted GZMB from
NKG2D-synNotch pair when using the optimal values from (D) round 1, (E) round 11 and (F)
round 15 from (A – C). (G – I) The amount of secreted PRF1 from NKG2D-synNotch pair when
using the optimal values from (G) round 1, (H) round 11 and (I) round 15 from (A – C). (D – I)
Markers, mean model prediction; error bars, one standard deviation; lightest shade, negative
cooperativity (V = 0.5); intermediate shade, no cooperativity (V = 1); darkest shade, positive
cooperativity(V=1.5).
123
FigureS3-15|Thedepletionofphospho-CD16limitsthemagnitudeofsteadystateGZMB
and PRF1. (A – D) The time series data of secreted GZMB when Rituximab is equal to (A) 10
-2
NM, (B) 10
-1
NM, (C) 1 NM and (D) 10 NM. (E – H) The time series data of secreted PRF1 when
Rituximab is equal to (E) 10
-2
NM, (F) 10
-1
NM, (G) 1 NM and (H) 10 NM. (I – L) The time series
dataofintracellularphospho-CD16whenRituximabisequalto(I)10
-2
NM,(J)10
-1
NM,(K)1NM
and (L) 10 NM. Solid line, mean model prediction; shaded area, one standard deviation; gray, no
control set (round 1 from Figure 4C); green, max control set (round 3 from Figure 4C); orange,
semi-maxcontrolset(round6fromFigure4C).
124
Figure S3-16 | The depletion of phospho-NKG2D limits the magnitude of steady state GZMB and
PRF1.(A–D)ThetimeseriesdataofsecretedGZMBwhenMICAisequalto(A)10-2OM,(B)10-
1OM,(C)1OMand(D)10OM.(E–H)ThetimeseriesdataofsecretedPRF1whenMICAisequal
to(E)10-2OM,(F)10-1OM,(G)1OMand(H)10OM.(I–L)Thetimeseriesdataofintracellular
phospho-NKG2D when MICA is equal to (I) 10-2 OM, (J) 10-1 OM, (K) 1 OM and (L) 10 OM. Solid
line, mean model prediction; shaded area, one standard deviation; blue, no control set (round 1
fromFigure4F);green,maxcontrolset(round11fromFigure4F);orange,semi-maxcontrolset
(round15fromFigure4F).
125
Chapter 4: A population-balance model of natural killer cell and
tumor cell interactions
Portions of this chapter are adapted from a manuscript by:
Sahak Z. Makaryan and Stacey D. Finley. In preparation.
4.1 Abstract
The effectiveness of natural killer (NK) cell-based cancer immunotherapies relies on
eliminating cancer cells upon contact. In practice, however, such therapies can be limited in
efficacy since cancer cells can mutate and change their receptor expression profile to evade
immune detection. Excitingly, mathematical modeling can provide actionable insight into
strategies that mitigate immune evasion – especially when applying dynamic game-theoretic
concepts to NK and tumor cell dynamics. Here, we present a system of partial differential
equations (PDEs) to describe how the two cell types are changing in concentration and how
they are changing their receptor concentrations as a result of cell-to-cell interactions. The
system is solved by using the well-known finite volume method. Furthermore, we
demonstrate how the initial distribution of the NK cell population over the receptor
concentration-space and the tumor cell properties influence the dynamics of the tumor cell
population. In summary, we inform strategies that enhance the anti-tumor activity of NK
cells.
126
4.2 Introduction
Natural killer (NK) cells are important immune effector cells that protect the host
from diseased cells such as virally infected cells and cancer cells
1,2
. In fact, when NK cells
engage with diseased cells, which express stress-induced stimulatory antigens on their cell
surface, NK cells become activated and kill the diseased cell by secreting the cytolytic
molecules granzyme B (GZMB) and perforin-1 (PRF1)
27,77–79
. Specifically, PRF1 mediates the
formation of pores on the target cell membrane, enabling GZMB to infiltrate and induce
apoptosis. Furthermore, the secretion of cytolytic molecules is mediated by multiple NK cell
receptor signaling pathways
3
, including the natural cytotoxicity receptors (e.g., NKp46,
NKp30), 2B4 (CD244), DNAM-1 (CD226), CD16 and NKG2D. However, NK cells can be
inhibited from cell killing when their inhibitory receptors are engaged, such as the
CD94/NKG2A receptor and members of the killer-cell immunoglobulin-like receptors
(KIRs). Inhibitory receptors are crucial for proper NK cell function
142–144
as they prevent NK
cells from damaging normal, healthy tissue. In fact, normal cells express major
histocompatibility molecules of class I (MHC-I) and non-classical human leukocyte antigens
(HLAs), which act as inhibitory ligands to NK cells. It is hypothesized that inhibitory
receptors help tune the responsiveness of NK cells via the Rheostat model
145,146
by
precluding over-activation of NK cells and by lowering the threshold for activation. Ardolino
et al
17
showed that NK cells extracted from the tumor microenvironment of MHC-I-deficient
tumors in mice were hypo-responsive to ex vivo cytokine and antibody-mediated
stimulation. While counterintuitive, NK cell inhibitory receptors help enhance the
responsiveness of NK cells to stimulation.
127
NK cell receptors mediate activation or inhibition through adaptor molecules and
other signaling molecules. For most natural cytotoxicity receptors, as well as for the CD16
receptor, the adaptor molecule CD3𝜁 facilitates the activation of the downstream signaling
molecules by recruiting the signalosome – a supramolecular complex of proteins – to the
membrane when the receptors are engaged with antibodies, ligands or drugs
23,53,81
. This
sequence of reactions brings together the Src family of kinases (SFKs) to the membrane near
the receptor
39,86,147
. Then, SFKs phosphorylate, and thus activate, intermediary proteins and
enzymes that ultimately lead to cell activation. In actuality, SFKs are recruited the
immunoreceptor tyrosine-based activation motifs (ITAMs) of the intracellular domain of
CD3𝜁. The receptor NKG2D mediates cell activation through the DAP10 adaptor molecule.
Dissimilar to CD3𝜁, DAP10 does not have any ITAMs; instead, the molecule contains an
amino acid sequence capable of recruiting SFKs that is similar in structure to the co-
stimulatory molecule CD28 in T cells
32,45,82,83
. The inhibitory receptors CD94 and KIRs,
however, contain immunoreceptor tyrosine-based inhibitory motifs (ITIMs), which help
recruit the phosphatases SHP-1, SHP-2 and SHIP to the membrane upon receptor
stimulation
39,148–150
. Indeed, these phosphatases dephosphorylate, and thus prevent
activation, of the signalosome. According to the Rheostat model
145,146
, NK cells almost always
simultaneously experience both stimulatory and inhibitory signals in vivo, and NK cell
activation is determined by the balance between such signals. In a physiological context, NK
cells encounter stimulatory and inhibitory signals when engaging with tumor cells
94
. Indeed,
it has been demonstrated
151
that tumor cells employ a variety of mechanisms for escaping
immune detection
152,153
, including changes to their receptor concentration profiles on their
cell surfaces. By increasing the expression of inhibitory antigens (e.g., HLAs), or decreasing
128
the expression of stimulatory antigens (e.g., MICA), tumor cells can evade NK cell-mediated
detection and killing. In addition, tumor cells can co-opt the surrounding stromal cells to
produce and secrete inhibitory cytokines
94
that prevents NK cell activation. Moreover, such
cytokines lead to an increase in inhibitory receptor expression on the surface of NK cells and
a decrease in stimulatory receptor expression, which unfortunately has been correlated with
poor patient prognoses in cancer patients
93
.
Given the complexity underpinning tumor-immune cell interactions, researchers
have employed computational approaches to study the tumor-immune dynamics in order to
provide insight into therapeutic strategies
40,150,154
. For instance, researchers have developed
mathematical models NK
38,39,41,86,147
and T cell
66,155–157
signaling to attain insight into
enhancing immune cell activation for the purpose of tumor cell killing. Elsewhere, game
theoretic approaches
158–161
to tumor growth dynamics have gained the attention of scientists
and physicians alike due to their intuitive framework; that is, the tumor cell population is
always maximizing their total density while physicians are always trying to minimize the
total number of tumor cells via therapeutic agents. Indeed, this dynamic game theoretic
framework provides a natural setting to study immune evasion. However, the above studies
either model the intracellular environment of the immune cell or model the intercellular
interactions between immune cells and tumor cells, but neither both. This leaves the
question of how heterogeneous populations of immune cells and tumor cells – expressing
various concentrations of stimulatory and inhibitory receptors – interact with each other.
In this work, we present a mathematical model of NK and tumor cell dynamics,
comprised of a system of partial differential equations (PDEs) to describe how the densities
of each cell type changes over time as well as over their receptor concentration-space. The
129
system of PDEs is solved via the finite-volume method. We demonstrate how the initial
distribution of NK cells over their receptor concentration-space impacts the tumor cell
population. In total, we provide a quantitative framework that can be applied to better
understand how the receptor distribution of NK cells limits tumor growth.
4.3 Methods
4.3.1 Population-balance modeling
NK cells express a variety of stimulatory and inhibitory receptors that mediate cell
activation or inhibition, respectively, upon stimulation
3
. Moreover, since NK cell populations
are seldomly homogeneous in their receptor expression profiles
43
, we require a model to
quantify how the number of NK cells are evolving over time as well as how their receptor
concentrations are changing over time. We require this also for the tumor cell population.
Excitingly, population-balance models enable us to model such dynamics
162–164
by
representing the system with the well-known advection equation from the theory of partial
differential equations (PDEs). We supply the explicit model equations, including initial and
boundary conditions, in the sections below for each cell type separately and summarize the
general model dynamics below. Briefly, we represent the dynamics using a system of coupled
PDEs that take the following general form:
𝜕𝑧
𝜕𝑡
+∇∙(𝑋𝑧) =𝑓(𝑡,𝑥,𝑧).
130
The variable 𝑧 is an evolving density function defined on the receptor concentration-
space. In other words, 𝑧 =𝑧(𝑡,𝑥) represents the number of cells expressing a receptor with
concentration in the infinitesimal volume [𝑥,𝑥+d𝑥] and in the infinitesimal time interval
[𝑡,𝑡+d𝑡]. To make the receptor concentration-space explicit, we define Ω to be an 𝑁-cell;
that is, Ω= [0,1]×⋯×[0,1] (in units of 𝜇M) is the cartesian product of 𝑁 unit intervals, with
𝑁 signifying the (total) number of receptors expressed by NK cells, tumor cells as well as
their complexes. The second term on the left-hand side of the equation above is the
divergence operator acting on the (velocity) vector field 𝑋, scaled by the density function 𝑧.
The vector field 𝑋 dictates how the density function 𝑧 changes over the receptor
concentration-space at each time 𝑡. Specifically, 𝑋 =r
FT
7
F`
,…,
FT
8
F`
t, where 𝑥
5
is the 𝑖
th
receptor
Figure 4-1 | Model schematic of NK- and tumor-cell interactions. NK cells are divided into three
separate populations: active, idle and proliferating. Active NK cells can bind to tumor cells,
forming a cell complex. The complex can dissociate or can initiate signaling. If the stimulatory
signal is stronger than the inhibitory signal, NK cells will kill the tumor cell, and thus forming an
idle NK cell and a dead tumor cell. Black arrows, stimulation.
131
expressed by the cell. Here, we use mass-action kinetics to describe how the receptor
concentrations are changing with time. Hence, the vector field 𝑋, at each time 𝑡, attaches a
tangent vector at each point in the receptor concentration-space, pointing in the direction of
where the cell is moving. Roughly speaking, 𝑋 describes the flow of a system of ordinary
differential equations (ODEs) in receptor concentration-space. Here, we assume the physical
environment is homogeneous. The right-hand side of the equation above (𝑓) represents the
net source term for the density function, delineating how each cell is created and destroyed
(Figure 4-1).
4.3.2 NK cell population
We distinguish between active NK cells (𝑁𝐾) and idle NK cells (𝐼). We assume that
idle NK cells – which are produced after a NK cell kills a tumor cell – does not engage with
tumor cells but remains idle for some time until re-joining the active NK cell population
140
.
All subtypes of NK cells express the receptors NKG2D and CD94. The stimulatory receptor
NKG2D belongs to the family of C-type lectin-like receptors, and is known to bind to stress-
induced antigens expressed by cancer cells and other cells that have undergone
transformation
11,32,45
. The inhibitory receptor CD94 prevents cell activation when it binds to
the inhibitory ligands MHC-I and the non-classical human leukocyte antigens (HLA)
148,151,165
.
Thus, 𝑁𝐾 and 𝐼 are evolving density functions defined on a 4-cell; we denote this space as
Ω
[\
. It is easy to check that Ω
[\
is a subset of Ω by setting the remaining dimensions (i.e.,
receptors expressed by tumor cells and complexes) to zero, and as a result, we have Ω
[\
⊂
𝜕Ω (the boundary of Ω).
132
The dynamics of active (𝑁𝐾) NK cells are defined as:
𝜕𝑁𝐾
𝜕𝑡
+∇∙(𝑋
[\
𝑁𝐾) =−𝑘
FGH
𝑁𝐾−𝑘
'(
𝑁𝐾 𝑇 d𝐿
a
3
+𝑘
'00
𝐶 d𝐿
a
9
\a
:;
+𝑘
c,_$
𝐼.
Based on data from Srpan et al.
53
, the receptor concentrations of NK cells remain constant
when NK cells are left unstimulated. We interpret this to mean active NK cells change their
receptor concentrations only when interacting with tumor cells. Therefore, we set 𝑋
[\
=
(0,0) which effectively makes the equation above an ODE with respect to time. The first term
on the right-hand side of the above equation is the constant decay of NK cells
40,154
. The
second term is the loss of active NK cells due to complex formation with tumor cells, with a
constant rate (𝑘
'(
) taken from Hoffman et al
41
. Since each active NK cell is free to bind to any
and all tumor cells, we integrate the tumor cell density function (𝑇) over the entire domain
(Ω
!
) to arrive at the total number of tumor cells. Here, we specify d𝐿 to be the volume
measure or volume form on Ω
!
. Conversely, the third term on the right-hand side is the gain
of active NK cells by dissociation of the complex (𝐶) with constant rate 𝑘
'00
. The complex
density function 𝐶 is defined space Ω
&
=Ω
[\
⋃Ω
!
, and therefore inheriting the receptor
concentrations from both active NK cells and tumor cells alike. Moreover, we define 𝐶 to be
such that the cells are bound to each other, but the receptors are not (see section 4.3.4). The
fourth term is the gain of active NK cells from idle NK cells (𝐼), with constant rate 𝑘
c,_$
estimated by Choi et al
140
, after a killing event occurred.
The initial condition of active NK cells is given by 𝑁𝐾(0,𝑅
2
,𝑅
7
) =𝑁𝐾
1
(𝑅
2
,𝑅
7
), where
𝑅
2
=NKG2D, 𝑅
7
=CD94. We leave 𝑁𝐾
1
in this general form as we will investigate how the
133
initial distribution of active NK cells over the receptor concentration-space affects the tumor
cell dynamics. Furthermore, NK cells have been shown to possess variability in their killing
potential
77,139,140
, where each individual NK cell can kill only a finite number of targets.
Indeed, this feature of NK cells is yet another source of heterogeneity; however, unlike the
receptor concentrations, this killing potential is discrete in nature. We model the killing
potential by fractioning the NK cell population into discrete sub-populations based on data
collected by Choi et al
140
. As an example, a given NK cell having a killing potential of 𝑖 will
move to the 𝑖−1 sub-population once it kills a single tumor cell. The dynamics and initial
condition of 𝑁𝐾 as defined above are left unchanged; we simply add a subscript 𝑖 to denote
the sub-population.
The dynamics of the idle NK cells are given by:
𝜕𝐼
𝜕𝑡
+∇∙(𝑋
d
𝐼)=−𝑘
FG+,e
𝐼−𝑘
c,_$
𝐼+𝑘
f_**
𝐻(𝐸
2
−𝐸
7
)𝐾 d𝐿d𝐸
a
;
\a
:;
.
Similar to the active NK cells, we assume idle NK cells are not changing their receptor
concentrations during their period of idleness
53,140
, implying 𝑋
d
= (0,0) and making the
above equation an ODE with respect to time. The first term on the right-hand side reflects
the death of idle NK cells after tumor cell killing
53
; this is referred to as activation-induced
death. The second term denotes the rate at which idle NK cells join the active NK cell
population
140
. The final term on the right-hand side is the production of idle NK cells after a
killing event has occurred. The quantity 𝐾 is the signaling complex formed by both active NK
cells and tumor cells where the receptors are engaged and actively signaling – as opposed to
134
the complex 𝐶 where the cells are merely bound to each other, but the receptors are not
engaged (see section 4.3.4). The quantity 𝐻 is the Heaviside function, with 𝐸
2
=𝑅
2
:𝐿
2
(NKG2D complexed with MICA) and 𝐸
7
=𝑅
7
:𝐿
7
(CD94 complexed with HLA). Note that 𝐸
2
is stimulatory whereas 𝐸
7
is inhibitory. If 𝐸
2
>𝐸
7
, then 𝐻 =1 and the NK cell will kill the
tumor cell with constant rate 𝑘
f_**
(estimated by Choi, et al
140
); otherwise, 𝐻 =0 since there
is more inhibitory signal than stimulatory signal. Lastly, we have the initial condition of
𝐼(0,𝑅
2
,𝑅
7
) =0.
4.3.3 Tumor cell population
Similar to the NK cell populations, the tumor cell population (𝑇) is heterogeneous
with regards to their receptor expression profile. In this work, we assume tumor cells
express MICA (a stimulatory ligand that binds to the NK cell receptor NKG2D) and HLA (an
inhibitory ligand that binds to CD94 on NK cells). The dynamics of the tumor cell population
is given below:
𝜕𝑇
𝜕𝑡
+∇∙(𝑋
!
𝑇)=−ℎ𝑇+2ℎ∙ℙ(𝐿)∙
∫ 𝑇 d𝐿
a
3
Vol(Ω
!
)
−𝑘
'(
𝑇 𝑁𝐾 d𝑅
a
:;
+𝑘
'00
𝐶 d𝑅
a
9
\a
3
.
Here, 𝑋
!
= r
F=
6
F`
,0t, where 𝐿
2
=MICA and 𝐿
7
=HLA. The second component in the
vector field 𝑋
!
is zero since we assume tumor cells are not changing their expression of the
inhibitory ligand HLA protein. However, there is evidence to suggest tumor cells can shed
135
their expression of the stimulatory ligand MICA
152,166,167
. This is reflected in the receptor
dynamics :
F=
6<
F`
=−𝛼
g%GF
𝐿
2
. The 𝛼
g%GF
parameter was taken from Pak et al
166
.
The first term on the right-hand side of the dynamics represents the loss of tumor
cells due to proliferation. The function ℎ =𝑘
H)'c
r1−𝑘
+,"
∫ 𝑇 d𝐿
a
3
t is a logistic growth rate
that is dependent on the total number of tumor cells at time 𝑡. The parameters 𝑘
H)'c
and 𝑘
+,"
denote the constant growth rate (in units of min
-1
) and the reciprocal carrying capacity (in
units of cells
-1
), which we took from de Pillis et al
154
. The second term represents the gain of
new daughter cells due to proliferation. Similarly, the function ℎ is the logistic growth rate
as defined above with the constant 2 denoting binary fission and the ratio on right is the total
tumor cell density. The quantity ℙ(𝐿) is a (joint) probability measure defined on Ω
!
,
denoting what fraction of the tumor cell population that can be produce daughter cells with
receptor concentrations less than or equal to 𝐿. For simplicity, we assume the probability of
the coordinates 𝐿
2
and 𝐿
7
are independent. Here, we study two distinct cases where tumor
cells are either optimizing for (1) growth versus (2) evasion. For the former scenario, we set
ℙ(𝑥) =∫ 𝛿(𝑦) d𝑦
T
1
, for each tumor cell receptor, as this case allows tumor cells to generate
offspring with any possible combination of receptor concentrations. Alternatively, for the
evasion scenario, we set ℙ(𝐿
2
) =expx−(𝜆𝐿
2
)y and ℙ(𝐿
7
) =expx𝜆(𝐿
7
−1)y with 𝜆 =10;
this construction of the probability measure on Ω
!
forces a drift towards a high expression
of the inhibitory ligand HLA and a low expression of the stimulatory ligand MICA. The third
and fourth terms on the right-hand side are defined similarly in the case of active NK cells;
that is, the loss of tumor cells due to association with active NK cells and the gain of tumor
cells from dissociation of the non-signaling complex.
136
Since the continuous uniform distribution maximizes the uncertainty of an unknown
probability distribution, we assume the tumor cell population is initially uniformly
distributed over their respective receptor concentration-space (Ω
!
). This assumption
maximizes the heterogeneity of the tumor cell population. Thus, we have the following for
the initial distribution: 𝑇
1
(𝐿)=1 𝑁
1
⁄ for all 𝐿 ∈Ω
!
and 𝑁
1
is the initial number of tumor
cells. As for the boundary condition, we impose a no flux condition: −(𝑋
!
∙n)𝑇 =0 for all 𝐿 ∈
∂Ω
!
and n is the (outward) unit normal vector field.
4.3.4 Complex population
The complexes formed by active NK cells and tumor cells inherit the receptor
concentrations from both cell types. Furthermore, we distinguish between complexes that
are signaling (𝐾) from those that are not (𝐶). More specifically, signaling complexes (𝐾) have
at least one receptor complex formed (i.e., 𝐸
5
>0 for some 𝑖 ∈{1,2}) while the non-signaling
complexes (𝐶) do not have any receptor complex formed (i.e., 𝐸
5
=0 for all 𝑖 ∈ {1,2}). This
assumption is based on imaging studies
60,77,82
, showing NK cells binding to tumor cells before
cell lysis occurs. By definition, we have Ω
&
=Ω
[\
⋃Ω
!
and Ω
&
⊂Ω
\
with Ω
\
⊂𝜕Ω.
The dynamics of the non-signaling complexes (𝐶) are given below:
𝜕𝐶
𝜕𝑡
+∇∙(𝑋
&
𝐶) =𝑘
'(
𝑆𝑇−𝑘
'00
𝐶.
By definition, 𝑋
&
≡0 and 𝐶
1
(𝑅,𝐿) =0 for all (𝑅,𝐿) ∈Ω
&
initially. The dynamics of the
signaling complexes (𝐾) is defined as:
137
𝜕𝐾
𝜕𝑡
+∇∙(𝑋
\
𝐾) =−𝑘
f_**
𝐻(𝐸
2
−𝐸
7
)𝐾.
The right-hand side of the above equation has been described previously (see section 4.3.3).
Observe that the signaling cell complexes (𝐾) do not dissociate unless all receptor complexes
dissociate (𝐸
5
=0 for all 𝑖), meaning the signaling complex must first become non-signaling
to then dissociate. Here, the vector field 𝑋
\
= r
FW
6
F`
,
FW
=
F`
,
F=
6
F`
,
F=
=
F`
,
Fh
6
F`
,
Fh
=
F`
t. The separate
receptor axes are defined as follows. The NKG2D and MICA axis is given by
d𝑅
2
d𝑡
=𝛼
'00
𝐸
2
−𝛼
'(
𝑅
2
𝐿
2
=
d𝐿
2
d𝑡
=−
d𝐸
2
d𝑡
;
and the CD94 and HLA axis is given by
d𝑅
7
d𝑡
=𝛼
'00
𝐸
7
−𝛼
'(
𝑅
7
𝐿
7
=
d𝐿
7
d𝑡
=−
d𝐸
7
d𝑡
.
The parameters for each receptor axis can be found in Table 4-2. As far as the initial
condition, 𝐾
1
(𝑅,𝐿,𝐸)=0 for all (𝑅,𝐿,𝐸) ∈Ω
\
. We impose a mixed boundary condition.
Recall that Ω
&
is a subset of the boundary of Ω
\
by taking the 𝐸
5
’s to be zero. Hence, Ω
&
⊂
𝜕Ω
\
, and we define the Dirichlet boundary condition 𝐾(𝑡,𝑅,𝐿,0) =𝐶(𝑡,𝑅,𝐿) for all (𝑅,𝐿) ∈
Ω
&
and for all 𝑡 >0. This condition forces NK and tumor cells to bind together before any
138
receptor signaling occurs. Also, for the boundary 𝜕Ω
\
\ Ω
&
, we have the no flux condition:
−(𝑋
\
∙n)𝐾 =0 for all 𝑡 >0.
Table 4-1. List of cellular parameters
Parameter Value Units Description
𝑘
'(
3.56×10
J?
(cells×min)
J2
Association rate
41
𝑘
'00
5.97×10
J7
min
J2
Dissociation rate
41
𝑘
FGH
2.86×10
J<
min
J2
Degradation rate of active NK cells
40
𝑘
FG+,e
5.21×10
JK
min
J2
Activation-induced death rate
53
𝑘
f_**
1.00×10
J2
min
J2
Rate of tumor cell killing
140
𝑘
c,_$
1.00×10
J2
min
J2
Rate of idleness of NK cells
140
𝑘
+,"
2.33×10
J3
cells
J2
Reciprocal of carrying capacity
40
𝑘
H)'c
4.04×10
JK
min
J2
Growth rate of tumor cells
40
Table 4-2. List of molecular parameters
NKG2D-MICA axis Value Units Description
𝛼
'(
4.05×10
1
(µM×min)
J2
Association rate
168
𝛼
'00
2.34×10
1
min
J2
Dissociation rate
168
𝛼
g%GF
6.67×10
J4
min
J2
Shedding rate
166
CD94-HLA axis Value Units Description
𝛼
'(
1.08×10
2
(µM×min)
J2
Association rate
169
𝛼
'00
2.52×10
2
min
J2
Dissociation rate
169
139
4.4 Results
4.4.1 The initial rate of tumor growth is slower if the tumor selects for evasion
We constructed a system of PDEs to describe how NK cells and tumor cells evolve
over time and over their receptor concentration space. Since we are interested in how the
initial receptor distribution of NK cells impacts tumor growth, we solved the system of PDEs
with two different initial distributions of the NK cell population: (1) NK cells are normally
distribution over the NKG2D- and CD94-space and (2) NK cells express a high concentration
of the stimulatory receptor NKG2D and a low concentration of inhibitory receptor CD94
(Figure 4-2). In addition, we investigated in silico the cases where the tumor cells are either
optimizing for growth or evasion of NK cells (see section 4.3.3 and Figure 4-3). Taken
Figure 4-2 | Initial distribution of NK cells. (A) NK cells are normally distributed across the
NKG2D- and CD94-space. (B) NK cells are distribution of the receptor-space such that they
express a high level of NKG2D and a low level of CD94.
140
together, we solve the PDE system using the finite-volume method for these four cases
separately where both populations have one thousand cells initially. Surprisingly, the model
predicts that the initial receptor distribution of NK cells is less significant in inhibiting tumor
growth – at least in the short-term (Figure 4-4). Indeed, we observed that the growth of the
tumor is less dependent on the initial receptor distribution of the NK cells (compare Figure
4-4A to Figure 4-4B) but is more dependent on whether the tumor cells are either optimizing
for growth or evasion (compare Figure 4-4A to Figure 4-4C). Overall, we found that the initial
dynamics of the tumor cell population is more dependent on whether tumor cells are trying
to evade NK cell detection or trying to maximize their total density.
Figure 4-3 | Probability measure on 𝛀
𝐓
. The probability landscape of tumor cell proliferation. (A)
The tumor cell population can generate offspring with any possible combination of MICA and HLA
while in (B) the tumor cell population produces offspring with a high level of HLA expression and
a low level of MICA expression.
141
4.5 Discussion
A mathematical model of NK cell and tumor cell dynamics was presented, describing
how the populations change in their concentrations in addition to changing their receptor
concentration profiles. Indeed, the model consists of a system of PDEs, with the spatial
domain representing the receptor concentration-space. The molecular and cellular
parameters were taken from literature (see Table 4-1 and Table 4-2). We approximated the
solution to the system using the finite-volume approach. Furthermore, the model was
simulated to determine how initial receptor distribution of NK cells impacts tumor growth.
Figure 4-4 | Initial dynamics of NK and tumor cells. (A) NK cells are normally distributed and
tumor cells are optimizing for growth. (B) NK cells are distributed according to Figure 4-2B and
tumor cells are optimizing for growth. (C) NK cells are normally distributed and tumor cells are
optimizing for evasion. (D) NK cells are distributed according to Figure 4-2B and tumor cells are
optimizing for evasion.
142
Interestingly, the model predicts that the initial receptor distribution of NK cells is not a
significant determinant of the anti-tumor response mediated by NK cells; instead, we found
that the objective of the tumor cells (i.e., maximizing growth versus evasion) is more
influential in determining the rate of tumor growth. These simulations seem to imply the
effectiveness of NK cell-based therapies – at least in the initial dynamics – are dependent on
whether the tumor cells are trying to evade or grow.
The presented model is not without limitations, however. For example, the model
does not include an exhaustive list of receptors expressed by both NK cells and tumor cells
alike. Certainly, future work can address this limitation by expanding the molecular
dynamics of NK cells to include the natural cytotoxicity receptors (e.g., NKp30) as well as
their downstream signaling molecules. Such additions can represent a biologically more
realistic description of NK-tumor cell dynamics. Moreover, the influence of cytokines on NK
and tumor cells was not described in our current work. Also, we only analyzed the initial
dynamics; researchers in the future can adopt the model and simulate for a longer time
course. Fortunately, such limitations can be addressed by modifying the model accordingly.
Limitations notwithstanding, the current work provides a quantitative framework for
understanding how the NK cell and tumor cell populations – as well as how their receptor
concentrations – evolve over time. By following the works of Ramkrishna
162–164
, the
population-balance modeling method is apt approach for describing the evolution of such a
system. Indeed, we have applied these concepts to our model of NK cell and tumor cell
dynamics and demonstrated which factors accelerate tumor growth. In summary, we
provided a computational framework for generating actionable insight into enhancing the
anti-tumor efficacy of NK cell-based therapies.
143
4.6 Acknowledgements
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).
144
Chapter 5: Conclusions
5.1 Summary
The computational work presented here provided new mechanistic insight into
processes by which NK cells can be modified to promote and enhance activation. We first
developed a mathematical model that describes the NK cell signaling network mediated by
three, well-studied stimulatory receptors: CD16, NKG2D and 2B4. These three receptors
have different signaling motifs that mediate cell activation, with various kinetics. For
example, both CD16 and NKG2D yield more phosphorylated (i.e., activated) signaling species
than 2B4 when all three receptors are stimulated with the same concentration of ligand. In
fact, the co-stimulation of CD16 and NKG2D required fewer ligands to achieve the same level
of species activation as the co-stimulation of all three receptors. These findings indicate the
NK cell signaling network is more responsive to CD16 and NKG2D stimulation whereas the
stimulation of 2B4 produces a weaker signal. Furthermore, we found that the availability of
ligand determines which strategy leads to more activation of the signaling molecules.
Specifically, when the ligand concentration is low, it is better to engineer proteolytic-
resistant receptors in order to continue signaling to the downstream molecules. Conversely,
when the amount of ligand is abundant in the system, the inhibition of the phosphatase
unbridles signal transduction, which enables the Src family of kinases to activate the
signalosome and the downstream molecular species.
In addition to enhancing the activation of the signaling network, we next investigated
in silico how the NK cell signaling pathways can be modified to enable the continued
145
secretion of the cytolytic molecules granzyme B (GZMB) and perforin-1 (PRF1). The
motivation behind this work was to predict mechanisms that counteract the effects of NK
cell exhaustion – a common phenotype observed when NK cells experience repeated
receptor stimulation. To that effect, we expanded the NK cell signaling model to include the
dynamics of the cytolytic molecules. After fitting the model predictions to experimental data,
we performed a sensitivity analysis and found that phosphatase inhibition increases the
amount of secreted cytolytic molecules; however, this came with a cost of depleting the
intracellular pool of said molecules. Inspired by modern synthetic biological methods, we
appended the model to include the dynamics of the synthetic Notch (synNotch) signaling
pathway to promote the production, and to expedite the secretion, of the cytolytic molecules.
By applying concepts from optimal control theory to the NK cell degranulation model, we
observed that the optimal amount of synthetic material depends on the frequency of
receptor stimulation. Since it has been observed that NK cells can kill multiple target cells,
our stimulations suggest it is best to incorporate a plasmid vector to promote the production
of GZMB and PRF1 in addition to a high affinity DNA binding domain to accelerate
transcription.
In conjunction with the intracellular modification of NK cells to enhance the secretion
of their cytolytic molecules, we inquired into the dynamics of NK and tumor cell populations.
Specifically, we constructed a mathematical model to delineate how individual NK cells and
tumor cells interact and how each cell type alters their surface receptors as a consequence
of such interactions. Here, applied a population-balance modeling approach to describe how
NK cells and tumor cells evolve over time and how their receptors change in time due to cell
contact. In effect, a system of partial differential equations (PDEs) was developed to address
146
this question, where the spatial domain represented the receptor concentration-space. We
the finite-volume method to approximate the solution of the system of PDEs. Indeed, we
simulated the model and found that the initial receptor distribution of NK cells less impactful
on constraining tumor growth. Interestingly, we observed that the rate of tumor growth is
faster if the tumor cell population is maximizing for growth as opposed to evasion.
5.2 Future directions
The series of models presented here are useful in providing actionable insight for
optimizing NK cell-based immunotherapies; however, they are not fully representative of the
true complexity underpinning both the molecular and cellular dynamics of NK cells. Surely
the intracellular signaling model can incorporate more receptor pathways than CD16, 2B4
and NKG2D. It would be of interest to both systems biologists and experimentalists on how
NK cell activation is influenced by the presence of both stimulatory and inhibitory signals.
This insight could lead to new strategies for engineering NK cells to overcome the inhibitory
signals conveyed by tumor cells or by the tumor microenvironment and thereby mitigate
immune evasion or immunosuppression. Indeed, the aforementioned model can be adapted
to address such questions by including the dynamics of such inhibitory receptors. Then, it
may be possible to define a quantitative relationship between the stimulatory and inhibitory
receptors in regard to cell activation.
The question of mitigating NK cell exhaustion is of particular interest to biomedical
scientists and researchers, especially those interested in maximizing the effectiveness of NK
cell-based therapies. Although NK cell exhaustion appears to be mediated by several factors,
147
it is possible to engineer them to become robust to receptor stimulation. Indeed, our in silico
studies and simulations suggest NK cell exhaustion can be indirectly mitigated by including
a synthetic pathway such that NK cells replenish their intracellular reservoir of cytolytic
molecules when their stimulatory receptors become activated – in other words, when a
danger signal is present in the milieu. An experimental implementation of said pathway can
extend and validate this work, and therefore generate a robust population of NK cells capable
of killing multiple target cells. In addition, the synthetic pathway can include a downstream
plasmid to promote the production of the stimulatory receptors; this may prevent the
reduction in the density of stimulatory receptors on the cell surface of NK cells, which has
been correlated with poor prognoses in cancer patients. Overall, by replenishing the
intracellular pool of cytolytic molecules – and by replenishing the density of stimulatory
receptors – the NK cell may become robust to cell activation.
While the above strategies may enhance the performance of NK cells with respect to
cancer cell-killing, the molecular signaling network has a finite capacity for mediating cell
activation (i.e., the signaling molecules degrade over time). Thus, there is a need to promote
the proliferation and expansion of the NK cell population to restrict tumor growth. However,
given that proliferation and target cell-killing are mutually exclusive processes, it is not at all
obvious how to direct the NK cell population to minimize the mass of the tumor. The NK-
tumor cell model can be extended to address the aforementioned issue. Indeed, a dynamic
game-theoretic approach can be applied to the present model, where the tumor cells are
maximizing their total mass while the NK cells are minimizing the mass of the tumor. By
including extracellular molecules such as interleukin-15 (IL15) and monoclonal antibodies
– that tag tumor-associated antigens – to induce NK cell proliferation and target cell-killing,
148
respectively, the model can be interrogated to understand how much of each molecule can
is needed to minimize tumor growth. Moreover, this framework can be applied to a variety
of conditions to yield optimal strategies for different physiological scenarios such as the case
where NK cells are outnumbered by the tumor cells by (say) one order of magnitude, which
is often the case in late-stage cancer.
5.3 Conclusion
Ultimately, the work provides a framework with which to study NK cell activation in
addition to strategies that maximize their effectiveness with respect to cancer cell killing. It
is able to provide actionable insight as to how to engineer NK cells for optimal cell killing.
Additionally, it can he used to test and generate new hypotheses to guide future experiments.
This framework can be expanded upon in the future to make a complete model that can
predict the efficacy of NK cell-based immunotherapies and contributing a more quantitative
understanding of NK cell biology.
149
References:
1. Bryceson, Y. T. & Ljunggren, H.-G. Natural Killer Cells: Biology, Physiology and
Medicine – Part 1. J. Innate Immun. 3, 213–215 (2011).
2. Bryceson, Y. T. & Ljunggren, H.-G. Natural Killer Cells: Biology, Physiology and
Medicine – Part 2. J. Innate Immun. 3, 327–328 (2011).
3. Long, E. O., Sik Kim, H., Liu, D., Peterson, M. E. & Rajagopalan, S. Controlling Natural
Killer Cell Responses: Integration of Signals for Activation and Inhibition. Annu. Rev.
Immunol. 31, 227–258 (2013).
4. Nelson, M. H. & Paulos, C. M. Novel immunotherapies for hematological malignancies.
Immunol. Rev. 263, 90–105 (2015).
5. Dong, S. & Ghobrial, I. M. Immunotherapy for hematological malignancies. J. Life Sci.
Westlake Village Calif 1, 46–52 (2019).
6. Singh, A., Brito, I. & Lammerding, J. Beyond Tissue Stiffness and Bioadhesivity:
Advanced Biomaterials to Model Tumor Microenvironments and Drug Resistance.
Trends Cancer 4, 281–291 (2018).
7. Emon, B., Bauer, J., Jain, Y., Jung, B. & Saif, T. Biophysics of Tumor Microenvironment
and Cancer Metastasis - A Mini Review. Comput. Struct. Biotechnol. J. 16, 279–287
(2018).
8. Reid, S. E. et al. Tumor matrix stiffness promotes metastatic cancer cell interaction
with the endothelium. EMBO J. 36, 2373–2389 (2017).
9. Wullkopf, L. et al. Cancer cells’ ability to mechanically adjust to extracellular matrix
stiffness correlates with their invasive potential. Mol. Biol. Cell 29, 2378–2385 (2018).
150
10. Dong, M. & Blobe, G. C. Role of transforming growth factor-β in hematologic
malignancies. Blood 107, 4589–4596 (2006).
11. Castriconi, R. et al. Transforming growth factor β1 inhibits expression of NKp30 and
NKG2D receptors: Consequences for the NK-mediated killing of dendritic cells. Proc.
Natl. Acad. Sci. 100, 4120–4125 (2003).
12. Donatelli, S. S. et al. TGF-β–inducible microRNA-183 silences tumor-associated natural
killer cells. Proc. Natl. Acad. Sci. 111, 4203–4208 (2014).
13. Palazon, A., Goldrath, A. W., Nizet, V. & Johnson, R. S. HIF Transcription Factors,
Inflammation, and Immunity. Immunity 41, 518–528 (2014).
14. Balsamo, M. et al. Hypoxia downregulates the expression of activating receptors
involved in NK-cell-mediated target cell killing without affecting ADCC. Eur. J.
Immunol. 43, 2756–2764 (2013).
15. Sarkar, S. et al. Hypoxia Induced Impairment of NK Cell Cytotoxicity against Multiple
Myeloma Can Be Overcome by IL-2 Activation of the NK Cells. PLOS ONE 8, e64835
(2013).
16. Das, S., Varalakshmi, C. & Khar, A. Target-cell-induced anergy in natural killer cells:
suppression of cytotoxic function. Cancer Immunol. Immunother. 49, 109–115 (2000).
17. Ardolino, M. et al. Cytokine therapy reverses NK cell anergy in MHC-deficient tumors. J.
Clin. Invest. 124, 4781–4794 (2014).
18. Jewett, A. & Tseng, H.-C. Tumor Microenvironment may Shape the Function and
Phenotype of NK Cells Through the Induction of Split Anergy and Generation of
Regulatory NK Cells. in The Tumor Immunoenvironment (eds. Shurin, M. R., Umansky,
151
V. & Malyguine, A.) 361–381 (Springer Netherlands, 2013). doi:10.1007/978-94-007-
6217-6_15.
19. Bi, J. & Tian, Z. NK Cell Exhaustion. Front. Immunol. 8, (2017).
20. Sun, C. et al. High NKG2A expression contributes to NK cell exhaustion and predicts a
poor prognosis of patients with liver cancer. OncoImmunology 6, e1264562 (2017).
21. Kadri, N., Luu Thanh, T. & Höglund, P. Selection, tuning, and adaptation in mouse NK
cell education. Immunol. Rev. 267, 167–177 (2015).
22. Morvan, M. G. & Lanier, L. L. NK cells and cancer: you can teach innate cells new tricks.
Nat. Rev. Cancer 16, 7–19 (2016).
23. Wirthmueller, U., Kurosaki, T., Murakami, M. S. & Ravetch, J. V. Signal transduction by
Fc gamma RIII (CD16) is mediated through the gamma chain. J. Exp. Med. 175, 1381–
1390 (1992).
24. Long, E. O. Negative signaling by inhibitory receptors: the NK cell paradigm. Immunol.
Rev. 224, 70–84 (2008).
25. Rajasekaran, K. et al. Signaling in Effector Lymphocytes: Insights toward Safer
Immunotherapy. Front. Immunol. 7, (2016).
26. Kim, H. S., Das, A., Gross, C. C., Bryceson, Y. T. & Long, E. O. Synergistic Signals for
Natural Cytotoxicity Are Required to Overcome Inhibition by c-Cbl Ubiquitin Ligase.
Immunity 32, 175–186 (2010).
27. Bryceson, Y. T., March, M. E., Barber, D. F., Ljunggren, H.-G. & Long, E. O. Cytolytic
granule polarization and degranulation controlled by different receptors in resting NK
cells. J. Exp. Med. 202, 1001–1012 (2005).
152
28. Bryceson, Y. T., March, M. E., Ljunggren, H.-G. & Long, E. O. Synergy among receptors
on resting NK cells for the activation of natural cytotoxicity and cytokine secretion.
Blood 107, 159–166 (2006).
29. Smole, A., Lainšček, D., Bezeljak, U., Horvat, S. & Jerala, R. A Synthetic Mammalian
Therapeutic Gene Circuit for Sensing and Suppressing Inflammation. Mol. Ther. 25,
102–119 (2017).
30. Guillerey, C., Huntington, N. D. & Smyth, M. J. Targeting natural killer cells in cancer
immunotherapy. Nat. Immunol. 17, 1025–1036 (2016).
31. Theorell, J. et al. Immunomodulatory activity of commonly used drugs on Fc-receptor-
mediated human natural killer cell activation. Cancer Immunol. Immunother. 63, 627–
641 (2014).
32. Billadeau, D. D., Upshaw, J. L., Schoon, R. A., Dick, C. J. & Leibson, P. J. NKG2D-DAP10
triggers human NK cell–mediated killing via a Syk-independent regulatory pathway.
Nat. Immunol. 4, 557–564 (2003).
33. Ting, A. T., Karnitz, L. M. & Schoon, R. A. Fc7 Receptor Activation Induces the Tyrosine
Phosphorylation of Both Phospholipase C (PLC)-y1 and PLC-y2 in Natural Killer Cells.
5 (1992).
34. Kim, H. S. & Long, E. O. Complementary Phosphorylation Sites in the Adaptor Protein
SLP-76 Promote Synergistic Activation of Natural Killer Cells. Sci. Signal. 5, ra49–ra49
(2012).
35. Kwon, H.-J. et al. Stepwise phosphorylation of p65 promotes NF-κB activation and NK
cell responses during target cell recognition. Nat. Commun. 7, (2016).
153
36. Kwon, H.-J. & Kim, H. S. Signaling for Synergistic Activation of Natural Killer Cells.
Immune Netw. 12, 240 (2012).
37. Pappalardo, F., Palladini, A., Pennisi, M., Castiglione, F. & Motta, S. Mathematical and
Computational Models in Tumor Immunology. Math. Model. Nat. Phenom. 7, 186–203
(2012).
38. Das, J. Activation or Tolerance of Natural Killer Cells Is Modulated by Ligand Quality in
a Nonmonotonic Manner. Biophys. J. 99, 2028–2037 (2010).
39. Mesecke, S., Urlaub, D., Busch, H., Eils, R. & Watzl, C. Integration of Activating and
Inhibitory Receptor Signaling by Regulated Phosphorylation of Vav1 in Immune Cells.
Sci. Signal. 4, ra36–ra36 (2011).
40. Mahasa, K. J., Ouifki, R., Eladdadi, A. & Pillis, L. de. Mathematical model of tumor–
immune surveillance. J. Theor. Biol. 404, 312–330 (2016).
41. Hoffman, F. et al. A mathematical model of antibody-dependent cellular cytotoxicity
(ADCC). J. Theor. Biol. 436, 39–50 (2018).
42. Scherbakova, A., Lust, H., Everaus, H. & Aints, A. A mathematical model of natural killer
cell activity. Cytometry A 83A, 585–591 (2013).
43. Sternberg-Simon, M. & Mehr, R. MODELING NATURAL KILLER CELL REPERTOIRE
DEVELOPMENT AND ACTIVATION DYNAMICS. in BIOMAT 2012 132–147 (WORLD
SCIENTIFIC, 2013). doi:10.1142/9789814520829_0008.
44. Chen, R. et al. Molecular Dissection of 2B4 Signaling: Implications for Signal
Transduction by SLAM-Related Receptors. Mol. Cell. Biol. 24, 5144–5156 (2004).
154
45. Gilfillan, S., Ho, E. L., Cella, M., Yokoyama, W. M. & Colonna, M. NKG2D recruits two
distinct adapters to trigger NK cell activation and costimulation. Nat. Immunol. 3,
1150–1155 (2002).
46. Barda-Saad, M. et al. Cooperative interactions at the SLP-76 complex are critical for
actin polymerization. EMBO J. 29, 2315–2328 (2010).
47. Harris, L. A. et al. BioNetGen 2.2: advances in rule-based modeling. Bioinformatics 32,
3366–3368 (2016).
48. Szewczuk, L. M., Tarrant, M. K. & Cole, P. A. Protein Phosphorylation by Semisynthesis:
From Paper to Practice. Methods Enzymol. 462, 1–24 (2009).
49. Northrup, S. H. & Erickson, H. P. Kinetics of protein-protein association explained by
Brownian dynamics computer simulation. Proc. Natl. Acad. Sci. U. S. A. 89, 3338–3342
(1992).
50. Schlosshauer, M. & Baker, D. Realistic protein–protein association rates from a simple
diffusional model neglecting long-range interactions, free energy barriers, and
landscape ruggedness. Protein Sci. Publ. Protein Soc. 13, 1660–1669 (2004).
51. Kim, M.-S. et al. A draft map of the human proteome. Nature 509, 575–581 (2014).
52. Bar-Even, A. et al. The Moderately Efficient Enzyme: Evolutionary and
Physicochemical Trends Shaping Enzyme Parameters. Biochemistry 50, 4402–4410
(2011).
53. Srpan, K. et al. Shedding of CD16 disassembles the NK cell immune synapse and boosts
serial engagement of target cells. J. Cell Biol. 217, 3267–3283 (2018).
155
54. Dong, Z. et al. The Adaptor SAP Controls NK Cell Activation by Regulating the Enzymes
Vav-1 and SHIP-1 and by Enhancing Conjugates with Target Cells. Immunity 36, 974–
985 (2012).
55. Schneider, C. A., Rasband, W. S. & Eliceiri, K. W. NIH Image to ImageJ: 25 years of image
analysis. Nat. Methods 9, 671–675 (2012).
56. Stuart, A. M. Inverse problems: A Bayesian perspective. Acta Numer. 19, 451–559
(2010).
57. Ghasemi, O. et al. Bayesian parameter estimation for nonlinear modelling of biological
pathways. BMC Syst. Biol. 5, S9 (2011).
58. Binstadt, B. A. et al. SLP-76 Is a Direct Substrate of SHP-1 Recruited to Killer Cell
Inhibitory Receptors. J. Biol. Chem. 273, 27518–27523 (1998).
59. Pérez-Quintero, L.-A. et al. EAT-2, a SAP-like adaptor, controls NK cell activation
through phospholipase Cγ, Ca++, and Erk, leading to granule polarization. J. Exp. Med.
211, 727–742 (2014).
60. Barber, D. F., Faure, M. & Long, E. O. LFA-1 Contributes an Early Signal for NK Cell
Cytotoxicity. J. Immunol. 173, 3653–3659 (2004).
61. Rudin, W. Real and complex analysis. (McGraw-Hill, 1987).
62. Clarke, F. Functional analysis, calculus of variations and optimal control. (Springer,
2013).
63. Arthur, D. & Vassilvitskii, S. K-means++: the advantages of careful seeding. in In
Proceedings of the 18th Annual ACM-SIAM Symposium on Discrete Algorithms
(2007).
156
64. Pearson, K. LIII. On lines and planes of closest fit to systems of points in space. (1901)
doi:10.1080/14786440109462720.
65. Kwon, H.-J. et al. NK cell function triggered by multiple activating receptors is
negatively regulated by glycogen synthase kinase-3β. Cell. Signal. 27, 1731–1741
(2015).
66. Rohrs, J. A., Wang, P. & Finley, S. D. Predictive Model of Lymphocyte-Specific Protein
Tyrosine Kinase (LCK) Autoregulation. Cell. Mol. Bioeng. 9, 351–367 (2016).
67. Rohrs, J. A., Makaryan, S. Z. & Finley, S. D. Constructing Predictive Cancer Systems
Biology Models. bioRxiv 360800 (2018) doi:10.1101/360800.
68. Binstadt, B. A. et al. Sequential Involvement of Lck and SHP-1 with MHC-Recognizing
Receptors on NK Cells Inhibits FcR-Initiated Tyrosine Kinase Activation. Immunity 5,
629–638 (1996).
69. Frank, C. et al. Effective Dephosphorylation of Src Substrates by SHP-1. J. Biol. Chem.
279, 11375–11383 (2004).
70. Lorenz, U. SHP-1 and SHP-2 in T cells: two phosphatases functioning at many levels.
Immunol. Rev. 228, 342–359 (2009).
71. Altvater, B. et al. 2B4 (CD244) Signaling by Recombinant Antigen-specific Chimeric
Receptors Costimulates Natural Killer Cell Activation to Leukemia and Neuroblastoma
Cells. Clin. Cancer Res. 15, 4857–4866 (2009).
72. Li, Y., Hermanson, D. L., Moriarity, B. S. & Kaufman, D. S. Human iPSC-Derived Natural
Killer Cells Engineered with Chimeric Antigen Receptors Enhance Anti-tumor Activity.
Cell Stem Cell 23, 181-192.e5 (2018).
157
73. Glienke, W. et al. Advantages and applications of CAR-expressing natural killer cells.
Front. Pharmacol. 6, (2015).
74. Romanski, A. et al. CD19-CAR engineered NK-92 cells are sufficient to overcome NK
cell resistance in B-cell malignancies. J. Cell. Mol. Med. 20, 1287–1294 (2016).
75. Bollino, D. & Webb, T. J. Chimeric antigen receptor–engineered natural killer and
natural killer T cells for cancer immunotherapy. Transl. Res. 187, 32–43 (2017).
76. X, T. et al. Erratum: First-in-man clinical trial of CAR NK-92 cells: safety test of CD33-
CAR NK-92 cells in patients with relapsed and refractory acute myeloid leukemia. Am.
J. Cancer Res. 8, 1899–1899 (2018).
77. Prager, I. et al. NK cells switch from granzyme B to death receptor–mediated
cytotoxicity during serial killing. J. Exp. Med. jem.20181454 (2019)
doi:10.1084/jem.20181454.
78. Brown, A. C. N. et al. Remodelling of Cortical Actin Where Lytic Granules Dock at
Natural Killer Cell Immune Synapses Revealed by Super-Resolution Microscopy. PLoS
Biol. 9, e1001152 (2011).
79. Voskoboinik, I., Whisstock, J. C. & Trapani, J. A. Perforin and granzymes: function,
dysfunction and human pathology. Nat. Rev. Immunol. 15, 388–400 (2015).
80. Felices, M., Lenvik, T. R., Davis, Z. B., Miller, J. S. & Vallera, D. A. Generation of BiKEs and
TriKEs to Improve NK Cell-Mediated Targeting of Tumor Cells. in Natural Killer Cells
(ed. Somanchi, S. S.) vol. 1441 333–346 (Springer New York, 2016).
81. Schmohl, J. U., Gleason, M. K., Dougherty, P. R., Miller, J. S. & Vallera, D. A.
Heterodimeric Bispecific Single Chain Variable Fragments (scFv) Killer Engagers
158
(BiKEs) Enhance NK-cell Activity Against CD133+ Colorectal Cancer Cells. Target.
Oncol. 11, 353–361 (2016).
82. Bryceson, Y. T. & Ljunggren, H.-G. Tumor cell recognition by the NK cell activating
receptor NKG2D. Eur. J. Immunol. 38, 2957–2961 (2008).
83. Chang, Y.-H. et al. A Chimeric Receptor with NKG2D Specificity Enhances Natural Killer
Cell Activation and Killing of Tumor Cells. Cancer Res. 73, 1777–1786 (2013).
84. Benitez, A. C. et al. Expression, signaling proficiency, and stimulatory function of the
NKG2D lymphocyte receptor in human cancer cells. Proc. Natl. Acad. Sci. 108, 4081–
4086 (2011).
85. López-Larrea, C., Suárez-Alvarez, B., López-Soto, A., López-Vázquez, A. & Gonzalez, S.
The NKG2D receptor: sensing stressed cells. Trends Mol. Med. 14, 179–189 (2008).
86. Makaryan, S. Z. & Finley, S. D. Enhancing network activation in natural killer cells:
predictions from in silico modeling. 13.
87. Chen, X., Trivedi, P. P., Ge, B., Krzewski, K. & Strominger, J. L. Many NK cell receptors
activate ERK2 and JNK1 to trigger microtubule organizing center and granule
polarization and cytotoxicity. Proc. Natl. Acad. Sci. 104, 6329–6334 (2007).
88. Saborit-Villarroya, I. et al. The Adaptor Protein 3BP2 Binds Human CD244 and Links
this Receptor to Vav Signaling, ERK Activation, and NK Cell Killing. J. Immunol. 175,
4226–4235 (2005).
89. Song, J. J. et al. c-Cbl acts as a mediator of Src-induced activation of the PI3K-Akt signal
transduction pathway during TRAIL treatment. Cell. Signal. 22, 377–385 (2010).
90. Mace, E. M. et al. Cell biological steps and checkpoints in accessing NK cell cytotoxicity.
Immunol. Cell Biol. 92, 245–255 (2014).
159
91. Sanchez-Correa, B. et al. Decreased expression of DNAM-1 on NK cells from acute
myeloid leukemia patients. Immunol. Cell Biol. 90, 109–115 (2012).
92. Paul, S., Kulkarni, N., Shilpi & Lal, G. Intratumoral natural killer cells show reduced
effector and cytolytic properties and control the differentiation of effector Th1 cells.
OncoImmunology 5, e1235106 (2016).
93. Peng, Y.-P. et al. Comprehensive analysis of the percentage of surface receptors and
cytotoxic granules positive natural killer cells in patients with pancreatic cancer,
gastric cancer, and colorectal cancer. J. Transl. Med. 11, 262 (2013).
94. Makaryan, S. Z., Cess, C. G. & Finley, S. D. Modeling immune cell behavior across scales
in cancer. Wiley Interdiscip. Rev. Syst. Biol. Med. 12, e1484 (2020).
95. Bryceson, Y. T. & Long, E. O. Line of attack: NK cell specificity and integration of
signals. Curr. Opin. Immunol. 20, 344–352 (2008).
96. Geweke, J. & Tanizaki, H. Bayesian estimation of state-space models using the
Metropolis–Hastings algorithm within Gibbs sampling. Comput. Stat. Data Anal. 37,
151–170 (2001).
97. Mitra, E. D. & Hlavacek, W. S. Parameter estimation and uncertainty quantification for
systems biology models. Curr. Opin. Syst. Biol. 18, 9–18 (2019).
98. Lüdtke, N. et al. Information-theoretic sensitivity analysis: a general method for credit
assignment in complex networks. J. R. Soc. Interface 5, 223–235 (2008).
99. Kojima, R. & Fussenegger, M. Synthetic Biology: Engineering Mammalian Cells To
Control Cell-to-Cell Communication at Will. ChemBioChem 20, 994–1002 (2019).
100. P Teixeira, A. & Fussenegger, M. Engineering mammalian cells for disease diagnosis
and treatment. Curr. Opin. Biotechnol. 55, 87–94 (2019).
160
101. Xie, M. & Fussenegger, M. Designing cell function: assembly of synthetic gene circuits
for cell biology applications. Nat. Rev. Mol. Cell Biol. 19, 507–525 (2018).
102. Morsut, L. et al. Engineering Customized Cell Sensing and Response Behaviors Using
Synthetic Notch Receptors. Cell 164, 780–791 (2016).
103. Lin, C. & Yang, L. Long Noncoding RNA in Cancer: Wiring Signaling Circuitry. Trends
Cell Biol. 28, 287–301 (2018).
104. Morriss, G. R. & Cooper, T. A. Protein sequestration as a normal function of long
noncoding RNAs and a pathogenic mechanism of RNAs containing nucleotide repeat
expansions. Hum. Genet. 136, 1247–1263 (2017).
105. Salviano-Silva, A., Lobo-Alves, S., Almeida, R., Malheiros, D. & Petzl-Erler, M. Besides
Pathology: Long Non-Coding RNA in Cell and Tissue Homeostasis. Non-Coding RNA 4,
3 (2018).
106. Yao, R.-W., Wang, Y. & Chen, L.-L. Cellular functions of long noncoding RNAs. Nat. Cell
Biol. 21, 542–551 (2019).
107. Torres, M. et al. RNA Pull-down Procedure to Identify RNA Targets of a Long Non-
coding RNA. J. Vis. Exp. 57379 (2018) doi:10.3791/57379.
108. Zhang, X. et al. Mechanisms and Functions of Long Non-Coding RNAs at Multiple
Regulatory Levels. Int. J. Mol. Sci. 20, 5573 (2019).
109. Long Non-Coding RNAs: Methods and Protocols. vol. 1402 (Springer New York, 2016).
110. Shamir, M., Bar-On, Y., Phillips, R. & Milo, R. SnapShot: Timescales in Cell Biology. Cell
164, 1302-1302.e1 (2016).
161
111. Chan, L. Y., Mugler, C. F., Heinrich, S., Vallotton, P. & Weis, K. Non-invasive
measurement of mRNA decay reveals translation initiation as the major determinant
of mRNA stability. eLife 7, e32536 (2018).
112. Lugowski, A., Nicholson, B. & Rissland, O. S. Determining mRNA half-lives on a
transcriptome-wide scale. Methods 137, 90–98 (2018).
113. Audet, C. & Dennis, J. E. Mesh Adaptive Direct Search Algorithms for Constrained
Optimization. SIAM J. Optim. 17, 188–217 (2006).
114. Marsden, J. E., Sirovich, L. & John, F. Optimization: Algorithms and Consistent
Approximations. (Springer New York, 1997).
115. Optimization and control with applications. (Springer Science+Business Media, 2005).
116. Bertsimas, D., Gupta, V. & Kallus, N. Robust sample average approximation. Math.
Program. 171, 217–282 (2018).
117. Chen, X., Shapiro, A. & Sun, H. Convergence Analysis of Sample Average Approximation
of Two-Stage Stochastic Generalized Equations. SIAM J. Optim. 29, 135–161 (2019).
118. Kleywegt, A. J., Shapiro, A. & Homem-de-Mello, T. The Sample Average Approximation
Method for Stochastic Discrete Optimization. SIAM J. Optim. 12, 479–502 (2002).
119. Pagnoncelli, B. K., Ahmed, S. & Shapiro, A. Sample Average Approximation Method for
Chance Constrained Programming: Theory and Applications. J. Optim. Theory Appl.
142, 399–416 (2009).
120. Jain, N. K. et al. A high density CHO-S transient transfection system: Comparison of
ExpiCHO and Expi293. Protein Expr. Purif. 134, 38–46 (2017).
162
121. Schukur, L. & Fussenegger, M. Engineering of synthetic gene circuits for (re-)balancing
physiological processes in chronic diseases: Engineering of synthetic gene circuits.
Wiley Interdiscip. Rev. Syst. Biol. Med. 8, 402–422 (2016).
122. Cooper, J. A., Kaneko, T. & Li, S. S. C. Cell Regulation by Phosphotyrosine-Targeted
Ubiquitin Ligases. Mol. Cell. Biol. 35, 1886–1897 (2015).
123. Matalon, O. & Barda-Saad, M. Cbl ubiquitin ligases mediate the inhibition of natural
killer cell activity. Commun. Integr. Biol. 9, e1216739 (2016).
124. Ibarra, I. L. et al. Mechanistic insights into transcription factor cooperativity and its
impact on protein-phenotype interactions. Nat. Commun. 11, (2020).
125. Morgunova, E. & Taipale, J. Structural perspective of cooperative transcription factor
binding. Curr. Opin. Struct. Biol. 47, 1–8 (2017).
126. Park, J. et al. Dissecting the sharp response of a canonical developmental enhancer
reveals multiple sources of cooperativity. eLife 8,.
127. Estrada, J., Wong, F., DePace, A. & Gunawardena, J. Information Integration and Energy
Expenditure in Gene Regulation. Cell 166, 234–244 (2016).
128. Gregor, T., Tank, D. W., Wieschaus, E. F. & Bialek, W. Probing the Limits to Positional
Information. Cell 130, 153–164 (2007).
129. Das, J. Positive feedback produces broad distributions in maximum activation attained
within a narrow time window in stochastic biochemical reactions. J. Chem. Phys. 138,
015101 (2013).
130. Das, J., Kardar, M. & Chakraborty, A. K. Positive feedback regulation results in spatial
clustering and fast spreading of active signaling molecules on a cell membrane. J.
Chem. Phys. 130, 245102 (2009).
163
131. Lugagne, J.-B. et al. Balancing a genetic toggle switch by real-time feedback control and
periodic forcing. Nat. Commun. 8, (2017).
132. Vecchio, D. D. & Murray, R. M. Biomolecular Feedback Systems. 282.
133. Lodeiro, M. et al. The SHP-1 protein tyrosine phosphatase negatively modulates Akt
signaling in the ghrelin/GHSR1a system. Mol. Biol. Cell 22, 4182–4191 (2011).
134. Lu, W., Gong, D., Bar-Sagi, D. & Cole, P. A. Site-Specific Incorporation of a
Phosphotyrosine Mimetic Reveals a Role for Tyrosine Phosphorylation of SHP-2 in Cell
Signaling. Mol. Cell 8, 759–769 (2001).
135. Matalon, O. et al. Dephosphorylation of the adaptor LAT and phospholipase C–γ by
SHP-1 inhibits natural killer cell cytotoxicity. Sci. Signal. 9, ra54–ra54 (2016).
136. Belikov, S., Berg, O. G. & Wrange, Ö. Quantification of transcription factor-DNA binding
affinity in a living cell. Nucleic Acids Res. 44, 3045–3058 (2016).
137. Kribelbauer, J. F., Rastogi, C., Bussemaker, H. J. & Mann, R. S. Low-Affinity Binding Sites
and the Transcription Factor Specificity Paradox in Eukaryotes. Annu. Rev. Cell Dev.
Biol. 35, 357–379 (2019).
138. Rastogi, C. et al. Accurate and sensitive quantification of protein-DNA binding affinity.
Proc. Natl. Acad. Sci. 115, E3692–E3701 (2018).
139. Bhat, R. & Watzl, C. Serial Killing of Tumor Cells by Human Natural Killer Cells –
Enhancement by Therapeutic Antibodies. PLoS ONE 2, e326 (2007).
140. Choi, P. J. & Mitchison, T. J. Imaging burst kinetics and spatial coordination during
serial killing by single natural killer cells. Proc. Natl. Acad. Sci. 110, 6488–6493 (2013).
141. Romain, G. et al. Antibody Fc engineering improves frequency and promotes kinetic
boosting of serial killing mediated by NK cells. Blood 124, 3241–3249 (2014).
164
142. Shifrin, N., Raulet, D. H. & Ardolino, M. NK cell self tolerance, responsiveness and
missing self recognition. Semin. Immunol. 26, 138–144 (2014).
143. Huse, M., Catherine Milanoski, S. & Abeyweera, T. P. Building tolerance by dismantling
synapses: inhibitory receptor signaling in natural killer cells. Immunol. Rev. 251, 143–
153 (2013).
144. Jaeger, B. N. & Vivier, E. Natural Killer Cell Tolerance: Control by Self or Self-Control?
Cold Spring Harb. Perspect. Biol. 4, a007229–a007229 (2012).
145. Brodin, P., Kärre, K. & Höglund, P. NK cell education: not an on-off switch but a tunable
rheostat. Trends Immunol. 30, 143–149 (2009).
146. Joncker, N. T., Fernandez, N. C., Treiner, E., Vivier, E. & Raulet, D. H. NK Cell
Responsiveness Is Tuned Commensurate with the Number of Inhibitory Receptors for
Self-MHC Class I: The Rheostat Model. J. Immunol. 182, 4572–4580 (2009).
147. Makaryan, S. Z. & Finley, S. D. An optimal control approach for enhancing natural killer
cells’ secretion of cytolytic molecules. APL Bioeng. 4, 046107 (2020).
148. Almeida, C. R. et al. Human NK Cells Differ More in Their KIR2DL1-Dependent
Thresholds for HLA-Cw6-Mediated Inhibition than in Their Maximal Killing Capacity.
PLoS ONE 6, e24927 (2011).
149. Charoudeh, H. N. et al. Quantity of HLA-C surface expression and licensing of KIR2DL+
natural killer cells. Immunogenetics 64, 739–745 (2012).
150. Guillamón, C. F. et al. NK Cell Education in Tumor Immune Surveillance: DNAM-1/KIR
Receptor Ratios as Predictive Biomarkers for Solid Tumor Outcome. Cancer Immunol.
Res. 6, 1537–1547 (2018).
165
151. Campoli, M. & Ferrone, S. Tumor escape mechanisms: Potential role of soluble HLA
antigens and NK cells activating ligands. Tissue Antigens 72, 321–334 (2008).
152. Beatty, G. L. & Gladney, W. L. Immune Escape Mechanisms as a Guide for Cancer
Immunotherapy. Clin. Cancer Res. 21, 687–692 (2015).
153. Mohme, M., Riethdorf, S. & Pantel, K. Circulating and disseminated tumour cells —
mechanisms of immune surveillance and escape. Nat. Rev. Clin. Oncol. 14, 155–167
(2017).
154. de Pillis, L. G. & Radunskaya, A. E. Modeling Tumor–Immune Dynamics. in
Mathematical Models of Tumor-Immune System Dynamics (eds. Eladdadi, A., Kim, P. &
Mallet, D.) 59–108 (Springer, 2014). doi:10.1007/978-1-4939-1793-8_4.
155. Rohrs, J. A., Wang, P. & Finley, S. D. Understanding the Dynamics of T-Cell Activation in
Health and Disease Through the Lens of Computational Modeling. JCO Clin. Cancer
Inform. 3, 1–8 (2019).
156. Rohrs, J. A., Zheng, D., Graham, N. A., Wang, P. & Finley, S. D. Computational Model of
Chimeric Antigen Receptors Explains Site-Specific Phosphorylation Kinetics. Biophys.
J. 115, 1116–1129 (2018).
157. Rohrs, J. A., Siegler, E. L., Wang, P. & Finley, S. D. ERK Activation in CAR T Cells Is
Amplified by CD28-Mediated Increase in CD3ζ Phosphorylation. iScience 23, 101023
(2020).
158. Newton, P. K. & Ma, Y. Maximizing cooperation in the prisoner’s dilemma evolutionary
game via optimal control. bioRxiv 2020.07.13.201400 (2020)
doi:10.1101/2020.07.13.201400.
159. Staňková, K. Resistance games. Nat. Ecol. Evol. 3, 336–337 (2019).
166
160. Staňková, K., Brown, J. S., Dalton, W. S. & Gatenby, R. A. Optimizing Cancer Treatment
Using Game Theory: A Review. JAMA Oncol. 5, 96–103 (2019).
161. Cunningham, J. J., Brown, J. S., Gatenby, R. A. & Staňková, K. Optimal control to develop
therapeutic strategies for metastatic castrate resistant prostate cancer. J. Theor. Biol.
459, 67–78 (2018).
162. Ramkrishna, D. Population balances: theory and applications to particulate systems in
engineering. (Academic Press, 2000).
163. Ramkrishna, D. & Singh, M. R. Population Balance Modeling: Current Status and Future
Prospects. Annu. Rev. Chem. Biomol. Eng. 5, 123–146 (2014).
164. Shu, C.-C., Chatterjee, A., Dunny, G., Hu, W.-S. & Ramkrishna, D. Bistability versus
Bimodal Distributions in Gene Regulatory Processes from Population Balance. PLoS
Comput. Biol. 7, e1002140 (2011).
165. Allard, M. et al. Serum Soluble HLA-E in Melanoma: A New Potential Immune-Related
Marker in Cancer. PLoS ONE 6, (2011).
166. Pak, Y., Zhang, Y., Pastan, I. & Lee, B. Antigen shedding may improve efficiencies for
delivery of antibody-based anticancer agents in solid tumors. Cancer Res. 72, 3143–
3152 (2012).
167. Zhao, Y. et al. Prognostic value of MICA/B in cancers: a systematic review and meta-
analysis. Oncotarget 8, 96384–96395 (2017).
168. Li, P. et al. Complex structure of the activating immunoreceptor NKG2D and its MHC
class I–like ligand MICA. Nat. Immunol. 2, 443–451 (2001).
169. Kaiser, B. K., Pizarro, J. C., Kerns, J. & Strong, R. K. Structural basis for NKG2A/CD94
recognition of HLA-E. Proc. Natl. Acad. Sci. 105, 6696–6701 (2008).
167
Abstract (if available)
Abstract
The interactions between natural killer (NK) cells?a subset of effector immune cells?and tumor cells has gained traction over the past decade due to the innate tumor cell killing ability of NK cells. Indeed, researchers have developed many therapeutic strategies, ranging from antibody-directed methods (e.g., bi-specific engagers) to genetic modifications that promote NK cell activation upon contact with cancer cells. Unlike conventional chemotherapy, the efficacy of such cell-mediated therapies relies heavily on the NK cell’s ability to eliminate cancer cells in an effective and successful manner. There are a couple of non-negligible drawbacks with this approach: (1) NK cells are heterogenous both in their killing potential and in their expression of stimulatory receptors that mediate cancer cell killing and (2) the activation of NK cells through its stimulatory signaling network is nonlinear with respect to its input. Moreover, finding optimal therapeutic strategies requires a quantitative understanding of the aforementioned issues, and without such a framework our therapies will almost surely be sub-optimal. In this work, I address the above concerns by constructing and simulating a set of mathematical models that describes the inter- and intra-cellular signaling dynamics between NK cells and tumor cells. The models are applied to predict (1) which stimulatory receptors?alone or in combination?most potently activate NK cell signaling network, (2) which strategies enhance NK cell activation with respect to the availability of input, (3) which intracellular signaling species most strongly controls and regulates NK cell activation, (4) how to perturb the signaling network to unbridle cell activation, (5) how to modify NK cells?via synthetic biology approaches?such that they secrete more cytolytic molecules upon receptor stimulation and (6) whether the density of stimulatory receptors on the cell surface of NK cells is strong determinant of tumor growth inhibition. Ultimately, the application of such mathematical models can inform both clinically relevant therapeutic strategies as well as a deeper understanding of NK cell signaling and activation.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Mechanistic model of chimeric antigen receptor T cell activation
PDF
Improving antitumor efficacy of chimeric antigen receptor-engineered immune cell therapy with synthetic biology and combination therapy approaches
PDF
Mechanistic investigation of pro-angiogenic signaling and endothelial sprouting mediated by FGF and VEGF
PDF
Genome-scale modeling of macrophage activity in the colorectal cancer microenvironment
PDF
Microenvironmental and biomechanical regulation of mitochondrial membrane potential in cancer cells
PDF
Metabolomic and proteomic approaches to understanding senescence and aging in mammary epithelial cells
PDF
Understanding the impact of cell-to-cell variability on intracellular signaling in CAR cells through mathematical modeling
PDF
Mechanistic modeling of angiogenic factors network and cancer therapy
PDF
Understanding anti-angiogenic signaling and treatment for cancer through mechanistic modeling
PDF
Combination therapy for solid tumor
PDF
Membrane-bound regulation of hematopoietic stem cells
PDF
Tumor-immune agent-based modeling: drawing insights from learned spaces
PDF
Phospho-proteomic analysis of immune cell activation
PDF
Photoplethysmogram-based biomarker for assessing risk of vaso-occlusive crisis in sickle cell disease: machine learning approaches
PDF
System biology approaches to cancer and diabetes metabolism
PDF
Enhancing the specificity and cytotoxicity of chimeric antigen receptor Natural Killer cells
PDF
Roles of circadian clock genes in cancer
PDF
Computational investigation of cholinergic modulation of the hippocampus and directional encoding of touch in somatosensory cortex
PDF
Modeling anti-tumoral effects of drug-induced activation of the cell-extrinsic apoptotic pathway
PDF
Systems approaches to understanding metabolic vulnerabilities in cancer cells
Asset Metadata
Creator
Makaryan, Sahak Zack
(author)
Core Title
Insights from a computational approach to natural killer cell activation
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Biomedical Engineering
Degree Conferral Date
2021-08
Publication Date
07/26/2021
Defense Date
04/09/2021
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
Control,mathematical model,natural killer cells,OAI-PMH Harvest,optimization,Systems Biology
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Finley, Stacey (
committee chair
), Graham, Nicholas (
committee member
), Khoo, Michael (
committee member
), Shen, Keyue (
committee member
)
Creator Email
makaryan@usc.edu,makaryansahak@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-oUC15622959
Unique identifier
UC15622959
Legacy Identifier
etd-MakaryanSa-9886
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Makaryan, Sahak Zack
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Repository Email
cisadmin@lib.usc.edu
Tags
mathematical model
natural killer cells
optimization