# Modeling putative therapeutic implications of exosome exchange between tumor and immune cells

^{a}Center for Theoretical Biological Physics,- Departments of
^{d}Physics and Astronomy, ^{b}Chemistry, and^{e}Biosciences, Rice University, Houston, TX 77005-1827;^{c}Red and Charline McCombs Institute for the Early Detection and Treatment of Cancer, University of Texas MD Anderson Cancer Center, Houston, TX 77030; and^{f}School of Physics and Astronomy and The Sagol School of Neuroscience, Tel-Aviv University, Tel-Aviv 69978, Israel

See allHide authors and affiliations

Contributed by José N. Onuchic, August 29, 2014 (sent for review August 19, 2014; reviewed by Kenneth J. Pienta, Wolfgang Losert, and Sorin Solomon)

## Significance

A better understanding of mechanisms of immune evasion by cancer cells and the role of the tumor microenvironment is crucial for developing new effective cancer therapeutic strategies. The challenge is posed by the enormous complexity of both the immune system and the tumor microenvironment, and the intricate cancer–immunity signaling network. Here, we develop a tractable theoretical framework to study the putative role of exosome communication in the cancer–immunity interplay. Exosomes are small (30–200 nm) vesicles that transfer proteins, mRNAs, and microRNAs to nearby and faraway cells. Guided by this model, we compare the effectiveness of administering radiation therapy alone or in combination with immunotherapy, illustrating how the model can shed light on the design and assessments of combination therapies.

## Abstract

Development of effective strategies to mobilize the immune system as a therapeutic modality in cancer necessitates a better understanding of the contribution of the tumor microenvironment to the complex interplay between cancer cells and the immune response. Recently, effort has been directed at unraveling the functional role of exosomes and their cargo of messengers in this interplay. Exosomes are small vesicles (30–200 nm) that mediate local and long-range communication through the horizontal transfer of information, such as combinations of proteins, mRNAs and microRNAs. Here, we develop a tractable theoretical framework to study the putative role of exosome-mediated cell–cell communication in the cancer–immunity interplay. We reduce the complex interplay into a generic model whose three components are cancer cells, dendritic cells (consisting of precursor, immature, and mature types), and killer cells (consisting of cytotoxic T cells, helper T cells, effector B cells, and natural killer cells). The framework also incorporates the effects of exosome exchange on enhancement/reduction of cell maturation, proliferation, apoptosis, immune recognition, and activation/inhibition. We reveal tristability—possible existence of three cancer states: a low cancer load with intermediate immune level state, an intermediate cancer load with high immune level state, and a high cancer load with low immune-level state, and establish the corresponding effective landscape for the cancer–immunity network. We illustrate how the framework can contribute to the design and assessments of combination therapies.

Immunotherapeutic approaches have recently emerged as effective therapeutic modalities (1) exemplified by immune checkpoint blockade with anti–CTLA-4 to activate T-cells and induce tumor cell killing, which has been shown to be effective for some cancers but not others (2). A better understanding of the intricate interplay between cancer and the immune system, and of mechanisms of immune evasion and of hijacking of the host response by cancer cells, is relevant to the development of effective immunotherapeutic approaches (3⇓⇓–6).

The immune-based suppression of tumor development and progression is mediated through nonspecific innate immunity and antigen-specific adaptive immunity (7). However, cancer cells can inhibit the immune response, thus evading suppression in multiple ways (8) (see below for details), and additionally hijack the immune system to their advantage (3, 4). The challenge to understand the tumor–immune interplay stems from the dynamic nature, and complexity and heterogeneity of both the cancer cells and the immune system and their interactions through the tumor microenvironment (9).

Here we consider immune cells as consisting of macrophages (10), natural killer cells (11), cytotoxic T cells (12), helper T cells (13), and regulatory T cells (3). These various immune cells are produced, activated, and perform their functions separated by space and time, which contributes to additional complexity (14). Among the immune cells, dendritic cells (DCs) are the most efficient antigen-presenting cells to bridge innate immunity with adaptive immunity (15). DCs also secrete cytokines that promote the antitumor functions of both natural killer cells and macrophages (16, 17). We consider the tumor microenvironment as comprised of a heterogeneous population of cancer cells (18), stromal cells (19), and tumor-infiltrating immune cells (20). The interactions among these cell types contribute to tumor development and progression. Tumor-associated macrophages and cancer-associated fibroblasts regulate tumor metabolism and engender an immune-suppressive environment by secreting TGF-β and other cytokines (21). Fluctuations in energy sources and oxygen within a tumor contribute to malignant progression and cell phenotypic diversity (22, 23).

Though secreted factors play critical roles in cell–cell communications, here we focus on the additional role of cell–cell communication mediated by the exchange of special extracellular lipid vesicles called exosomes (24). These nanovesicles of ∼30–200 nm are formed in the multivesicular bodies and then released from the cell into the extracellular space (25). The exosomes carry a broad range of cargo, including proteins, microRNAs, mRNAs, and DNA fragments, to specific target cells at a remote location (26). Membrane markers assign the exosomes to specific targeted cells. Notably, upon entering the target cell, the exosomes induce modulation of cell function and even identity switch (phenotypic, epigenetic, and even genetic) (27). Exosomes have recently emerged as playing an important role in the immune system interaction with tumors (28, 29). Tumor-derived exosomes can promote metastatic niche formation by influencing bone marrow-derived cells toward a prometastatic phenotype through upregulation of c-Met (29). DCs have been shown to induce tumor cell killing through release of exosomes that contain potent tumor-suppressive factors such as TNF and through activation of natural killer cells, cytotoxic T cells, and helper T cells (24, 30⇓–32). However, tumor-derived exosomes (Tex) can directly inhibit the differentiation of DCs in bone marrow (33), which strongly inhibits the dendritic cell-mediated immune response to the tumor. In addition, Tex can also directly inhibit natural killer cells (34).

Mathematical models have been devised to study the complex interactions of cancer and immune system, including those that consider spatial heterogeneity (as reviewed in ref. 35) and those that consider spatially homogeneous populations (as reviewed in ref. 36). Cancer–immunity models have been constructed to investigate the effects of therapy (37⇓–39), cancer dormancy (40), and interactions with time delay (41). Other types of modeling methods have also been applied. For example, tumor growth has also been fitted to experimental data by artificial neural networks (42); a detailed network of cancer immune system has been modeled by multiple subset models (43).

In this study, we have developed an exosome exchange-based cancer–immunity interplay (ECI) model, to incorporate the special role of DCs and exosome-mediated communications. Distinct from the previous approaches, our modeling strategy is adapted from methodology used in studies of gene regulatory circuits, allowing us to check the multistability features of the system (44). We find that, by including exosome exchange, the cancer–immunity interplay can give rise to three quasi-stable cancer states, which may be associated with the elimination/equilibrium/escape phases proposed in the immunoediting theory (45). The ECI model is also capable of explaining tumorigenesis by considering the time evolution of immune responses. Guided by the treatment simulations, we assess the effectiveness of various therapeutic protocols with and without time delay and noise.

## Results

### Exosome Exchange-Based Cancer–Immunity Interplay Modeling Approach.

Here, we develop a minimal yet workable theoretical framework for modeling the cancer–immunity interplay. The concept is to devise ECI models that incorporate three generic (coarse-grained) components, i.e., the effective cancer cells (C), the dendritic cells (D), and the killer cells (K). C represents the cancer cells that actively interact with the immune system, e.g., those at the surface of a solid tumor. D represents precursor, immature, and mature DCs that present tumor-associated antigens. The cytotoxic T cells, effector B cells, helper T cells, and natural killer cells are lumped together and referred to as killer cells (K) because, as a group, they inhibit cancer cells after being activated by DCs.

To elucidate the effects of cell–cell communications among these representative cell types, the dynamics of C, D, and K cells are modeled by nonlinear ordinary differential equations. We consider two communication networks as illustrated in Fig. 1. For both cases, the C population is self-activated through exosome-mediated cell proliferation (46). The D population is also self-activated but by means of D-cell proliferation. D also activates proliferation of K (30⇓–32), and K induces maturation of D (47). Also, K targets C and induces C’s apoptosis (11⇓–13, 48, 49), and C inhibits K by slowing down its proliferation and increasing its apoptosis (50, 51). The two networks, however, differ from each other in the interactions between C and D. In case I (Fig. 1*A*), D directly kills C by exosomes (24) and/or some cytokines (52); C induces the D maturation by presenting tumor-associated antigens on the surface when the concentration of C is low (15, 53), but inhibits D differentiation by exosomes and/or cytokines when the concentration of C is above a threshold level (33, 54). The deterministic rate equations for the dynamics of C, D, and K densities are given by*SI Appendix*, Table S1. In case II (Fig. 1*B*), the exosome-mediated interactions are not included, thus the most prominent interaction between C and D is the activation of D by recognizing antigens that are associated with C. We remain using Eq. **1** to model the dynamics of the network, but change the parameters for the modified interactions, i.e., the fold changes

The effects of immune recognition are incorporated by varying the interactions from D to C and K to C using a recognition parameter (*ρ*) as follows:*ρ* represents the effective immune recognition: *ρ* = 0 corresponds to no recognition so that there is no effect from D or K to C, i.e., *ρ* = 1.0 corresponds to full recognition so that *ρ* can also be larger than 1, which represents the case of enhanced immune response above the basal level that can be induced, for example, by immunotherapy.

Immune recognition plays an essential role in tumor progression. We consider that nascent tumor cells are not initially recognized by the immune system at the very onset (*ρ* = 0). As the tumor continues to progresses, the *ρ* value gradually increases from 0 to 1 due to enhancement of its immune recognition. The evolving tumor can evade the immune system by multiple mechanisms, e.g., through mutations. During such a period, the *ρ* value gradually decreases. We discuss the effects of immune recognition on tumorigenesis. Notably, because the value of *ρ* increases at a slower pace than the dynamics of the cancer–immunity interplay, we study the latter for different constant level of *ρ* (quasi-static approximation). The framework may be further extended to include a dynamical equation for *ρ* that is coupled with the three dynamical equations for C, D, and K.

### Existence of Multiple Cancer States During Tumorigenesis.

We first consider two cases: case I in which exosome exchange is included, and case II in which exosome exchange is not included. For both cases, the quasi steady-state solutions of Eq. **1** are obtained by the phase-plane analysis (*Methods*).

In case I, when the immune recognition is low (*ρ* = 0.2), the interplay yields only one quasi-stable steady state, corresponding to a high cancer load with low immune level [denoted as the (H) or (1, 0) state; Fig. 2*A*]. The reason that this is the only possible state is that when the immune recognition is low, the effect of the immune system is not strong enough to limit cancer progression. In contrast, when immune recognition is high (*ρ* = 1.0), the interplay can give rise to three quasi-stable steady states; Fig. 2*B*), corresponding to the (H) state, an intermediate cancer load with high immune level [denoted as the (I) or (1/2, 1) state], and a low cancer load with intermediate immune level [denoted as the (L) or (0, 1/2) state]. More specifically, for *ρ* = 1.0, the effect of the immune system on cancer is comparable to the effect of cancer on the immune system, therefore allowing the existence of multiple quasi-stable states. The relative stability of each state can be inferred from the stochastic simulations (*Methods*). When *ρ* = 2.0, the immune system outweighs cancer, allowing only a low cancer load with low immune level (0, 0) state deterministically; it could, however, have some rare chances to switch to the high cancer load with low immune level state stochastically (*SI Appendix*, Fig. S1*A*). Notably, the predicted existence of three quasi-stable states in the ECI model is likely related to the clinical classification of the equilibrium stage during tumor progression as defined in the immunoediting theory (45), which considers both the host-protection and tumor-sculpting action of immune system on tumorigenesis.

In case II, the interplay does not include negative interactions between C and D. These interactions are mainly mediated by the exosome exchange. As a consequence, case II dynamics are significantly different from those of case I. When *ρ* = 0.2, the interplay leads to a high cancer load with high immune (1, 1) state (Fig. 2*C*). However, when *ρ* = 1.0, the interplay leads to bistability—existence of an intermediate cancer load with high immune level (1/2, 1) state and a low cancer load with intermediate immune level (0, 1/2) state, meaning that a third state does not exist (Fig. 2*D*). When *ρ* = 2.0, the interplay gives rise to a low cancer load with low immune level (0, 0) state in both the deterministic and the stochastic analyses (*SI Appendix*, Fig. S1*C*). In *SI Appendix*, we show a comparison with the stability analysis for two additional cases in which a part of exosome exchange interactions are not included (*SI Appendix*, Fig. S2).

The comparison between these different cases implies that exosome exchange facilitates tristability—the existence of the (H), (I), and (L) states. As we show next, the existence of the (I) state, which is possible due to the exosome exchange, has very important implications regarding the effectiveness of therapeutic strategies.

### Bifurcation Diagram for Immune Recognition.

To further assess the effect of immune recognition in case I, we calculate and investigate the bifurcation diagrams of the cancer–immunity interplay with respect to the level of *ρ* as a control parameter (Fig. 3). As is depicted in this figure, the bifurcation diagram exhibits four different regimes (phases), each corresponding to the existence of different quasi-stable cancer states: (*i*) a phase in which only the (L) state exists (denoted as {L}, blue); (*ii*) a phase in which both the (L) and (H) states exist (denoted as {L,H}, yellow); (*iii*) a phase in which the (L), (I), and (H) states exist (denoted as {L,I,H}, green); and (*iv*) a phase in which only the (H) state exists {denoted as {H}, red). Further inspection reveals that when *ρ* is above 1.1, the (L) state is the dominant one; when *ρ* is between 0.84 and 1.1, all three states exist, making the transitions among these states more likely; when *ρ* is below 0.84, the (H) state is the dominant one. It is worth noting that a bifurcation curve may cross itself, because the *y* axis of the bifurcation diagram is projected from a 3D phase-plane (Fig. 3 *C* and *D*); besides, the bifurcation curve for D (Fig. 3*B*) is more similar to that for the K (*SI Appendix*, Fig. S3) than that for C (Fig. 3*A*), showing that D correlates with K. For case II, the interplay allows the existence of, at most, two quasi-stable cancer states, and D correlates with C (*SI Appendix*, Fig. S4). The distinct dynamics of the two cases illustrate the importance of the exosome exchange.

For clarity, in the following, we focus on analysis for the network in case I.

### Immune Recognition and Tumorigenesis.

Before tumor onset, immune recognition has yet been established, i.e., *ρ* is zero, the C population is close to zero, and the D population is low. As the tumor starts to develop, it triggers immune recognition. The fate of tumor development is determined by the relative rate of tumor growth vs. the rate of increase in immune recognition. In general, for fast-growing tumors, the dynamics will lead to the (H) cancer state. More specifically, depending on how fast the immune system responds (the level of immune recognition increases), which is characterized by the mean lifetime (*ρ* = 1; *Methods*), the interplay (for the model parameters of cancer growth used here) can lead to the (H) state (when *A*. Further enhancement of immune recognition causes the transition of the tumor from the (H) or the (I) states to the (L) state (downward dashed arrows labeled as 1 and 2 in Fig. 3*A*), suggesting possible control of the tumor by the immune system and the putative importance of boosting the immune system during early stages of tumorigenesis.

### Effects of Tumor Evasion on Immune Recognition.

As is now widely recognized, cancer cells can evade immune suppression during tumor development, which causes the immune recognition to gradually decrease. In general, tumor size, as quantified by the size of the C population, increases when *ρ* decreases. Marked growth of C can be observed in some stages either deterministically (upward dashed arrows labeled as 3 and 4 in Fig. 3*A*), or stochastically when the system overcomes a dynamical barrier of the transition between two coexisting cancer states. Interestingly, such marked growth of C happens within phase {L,I,H} and the leftmost part of phase {L,H}. At these phases, the levels of both D and K are relatively higher than those at the other phases, which is consistent with experimental evidence (55). At the late tumor stage (low *ρ* again), certain processes, including epithelial–mesenchymal transition, may trigger cancer metastasis, in which some small clusters of cancer cells migrate to a remote location; from the new site, they initiate the development of a new tumor load from the metastatic cancer cells that may have different identity and therefore require renewed immune recognition for tumor control.

### Sensitivity to Model Parameters.

We further explore sensitivity of the existence of the intermediate state and the tristability (defined in *SI Appendix*, section 5) to model parameters assuming that each tumor is characterized by a specific set of parameters. For this aim, we calculated the bifurcation diagram with respect to *ρ* for different sets of parameters. More specifically, we randomly vary all 36 model parameters simultaneously away from the original values by a uniform distribution with maximum of range of ±*d* (chosen from 10%, 20% … 60%). We find that the intermediate (I) cancer state exists for a significant range of model parameters (blue columns in Fig. 4*A*). For example, the (I) state exists for ∼50% of the parameter sets (or 50% chance) when the 36 parameters are varied randomly within a range of 40% about their original values. The chance for observing the tristability (red columns in Fig. 4*A*) is slightly lower than that for the existence of the (I) state. Among those parameter sets in which the (I) state exists (blue columns in Fig. 4*B*), both the average range of existence for the (I) state (defined in *SI Appendix*, section 5) and its SD (*SI Appendix*, Fig. S6) go up with the width of the distribution from which the parameters are taken. Similar behavior is observed for the tristability (red columns in Fig. 4*B*). For a larger perturbation of the parameters, the range of existence for the (I) state is much wider than the range of existence for the tristability (Fig. 4*B*). The results imply that when exosome exchange is operative, the intermediate cancer state and the tristability are characteristic features of the ECI model for a wide range of model parameters.

We further assess the effect of individual parameters on the existence of the (I) state and the tristability (*SI Appendix*, Figs. S7 and S8). We find that when any one of the parameters is varied by 10%, the range of existence for the (I) state and of the tristability are only slightly affected (*SI Appendix*, Fig. S7 *A* and *B*), but a variation of 30% in a single parameter can sometimes lead to significant difference (*SI Appendix*, Fig. S7 *C* and *D*). The three most significant model parameters are associated with tumor growth (the Hill threshold of C self-activation *SI Appendix*, Fig. S9). These results imply that the existence of tristability and the (I) cancer state arises from a balance between cancer and immunity due to exosome exchange.

### Building the Cancer–Immunity Landscape.

Motivated by the notion of Waddington landscape for cell differentiation (56), we proceed to introduce an effective landscape (57) corresponding to the dynamical states of the cancer–immunity interplay; doing so helps to better understand the nature of the coexisting cancer states and to design and assess possible therapeutic protocols by visualization of their effect as trajectories connecting different cancer states (or as a transition rate problem) in established landscape. In principle, the dynamics should be presented as a 4D landscape whose axes are the densities of C, D, K, and the strength of the immune recognition *ρ*. Here, by assuming quasi-steady state for *ρ*, we consider 3D landscapes (C, D, and K) for different values of *ρ*. The landscape defines an effective potential computed as minus the logarithm of the probability [−log(P)] of each dynamical state (C, D, K) in the presence of noise. The noise could originate from multiple sources, including cell–cell communications, birth–death processes of individual cells, coarse-graining of different cell types, and the effect of some other cancer/immune cell types or factors that are not included in the ECI model.

In practice, it is easier to visualize the landscape in two dimensions. Therefore, we rely on projection of the 3D space onto a corresponding 2D one. For example, in Fig. 2 and in *SI Appendix*, Fig. S1, the landscape is projected onto a 2D space whose axes are C and D. As shown in *SI Appendix*, we find that the landscape could be better presented if we replace D with a combination of D and K as one of the axes (*SI Appendix*, Figs. S10 and S11), because D and K are highly correlated (although in a nonlinear manner). We propose that the concept of the cancer–immunity landscape is very valuable for the design and assessment of effectiveness of novel therapeutic strategies. In particular, the landscape makes it clear that the desired/optimal treatments are not those that simply reduce the tumor size. As illustrated further below, treatment effectiveness should be assessed by its efficacy to cross the dynamical barrier between high and low cancer states. It also shows that in most cases a direct transition from the (H) cancer state to the (L) state (H to L) is very hard to achieve. In the majority of cases, two-stage strategies are the desired ones—one treatment is used for effective H-to-I transition and a second treatment is used for effective I-to-L transition.

### Effects of Time Delay on the Cancer–Immunity Landscape.

Communication between C, D, and K cell types involves time delay. A longer time delay is expected for exosome exchange, because of the need for their production and transport across long distances. From a dynamical system perspective, time delay can have a significant effect on transitions between cancer states and the effect of noise on the cancer–immunity interplay. We note that steady-state solutions are not affected when the time delay is considered. However, the time delay might alter dynamics by changing the stability of the steady states (the cancer states) and the shape of the landscape. To this end, we evaluate the possible effects of time delays on the cancer–immunity interplay by performing simulations of the ECI model (case I)

### Evaluation of the Flow Toward Cancer–immunity States.

We examine the relative flow toward each of the stable states [(H), (I), and (L)] while incorporating the time delay in exosome exchange between C and D. From a dynamical systems perspective, we calculate, for each state, its corresponding basin of attraction—the initial conditions from which the cancer/immunity evolve toward the particular dynamical state. More specifically, assuming the subject is initially at the high cancer state (H), we simulate possible consequence of hypothetical treatments; this is done by changing the population values of both C and D, whereas K cells are not affected. Such a treatment changes the system from the (H) state to a new state in the phase plane (Fig. 5*A*). We then perform simulations with time delay (*Methods*) for the dynamics starting from the newly induced location (new values of C and D). The simulations reveal that the dynamics after the initial induced deviations generate a trajectory that ends in one of the quasi-stable states [(H), (I), and (L)]. We define the basin of attraction for the (H) state as the region in the landscape (the phase plane), from which the cancer–immunity evolves back toward the (H) state. The basins of attraction for the (I) and (L) states can be constructed in a similar manner.

### Effects of Time Delay on the Flow Toward Cancer–Immunity States.

Fig. 5*A*, *Upper Left*, shows the basins of attraction for the (H) (navy), (I) (blue), and (L) (light blue) states for case I and without time delay. An effective treatment should cause a transition from the (H) state to either the (I) or the (L) state. Our results indicate that if a treatment aims to only reduce the C population, it requires a reduction of more than 90%, otherwise the system is still within the basin of attraction of the (H) state. Alternatively, it is more effective if the treatment targets the D population or both the C and D populations. By adding and increasing time delay, the basin of attraction for the (H) state increases linearly for small values of time delay and is saturated for large values of time delay (Fig. 5*B*), and that for the (L) state is mostly unaffected.

### Model-Based Design and Assessment of Therapeutic Strategies.

We present a design of several hypothetical therapeutic strategies and assessment of their effectiveness. Again, we assume that a subject’s cancer is initially in the (H) state. Because the (H) state is located away from the (L) state in the phase space (Fig. 2*B*), it is hard to envision a treatment that can induce direct transitions from the (H) state to the (L) state (H2L). For this reason, in most cases (model parameter or cancer type) a more efficient protocol is likely to be based on a two-stage strategy: the first involves a treatment to induce the transition from the (H) state to the (I) state (H2I), and the second involves a treatment to induce the transition from the (I) state to the (L) state (I2L). We note that use of this more efficient two-stage strategy is suggested by the existence of the (I) state due to the exosome exchange.

To illustrate the concept in a quantitative manner, we mainly focus here on two types of common anticancer therapies: radiation therapy and DC-based immunotherapy. The effect of the radiation therapy is incorporated into the ECI model as elevation in apoptosis rates of C as well as of D and K during treatment. Note that apoptosis of D and K represents the inhibitory effect of radiation on the immune system (58, 59). DC immunotherapy is represented in the ECI model as elevation in proliferation rates of D cells (*Methods*).

### Assessment of Simulated Radiation Therapy.

The effect of hypothetical radiation therapy is shown in Fig. 6. This example is composed of 30 one-per-day radiation sessions from day 10 (point 0) to day 40 (point 1) (Fig. 6*B*). The population of D cells decreases during the 30 d (solid red line in Fig. 6*C*). After the therapy ends, the population of D cells increases and reaches maximum at approximately day 90 (point 2), and then gradually decreases toward its initial level before therapy started (point 3, dashed back line). These results are in excellent agreement with clinical data for breast cancer patients after 30-d radiation therapy (Fig. 6*D*; data from ref. 59). We note that the percentage of phagocytic activity of monocytes was measured in the clinical tests and not the density of D cells. However, the latter is correlated with the percentage of phagocytic activity of monocytes. We note that the model simulations are done independently of the clinical data, meaning that we do not fit the model parameters for this specific clinical test. For this reason, though in the model simulations the time recovery of the immune system to its value before the hypothetical radiation treatment is 7 wk, in the clinical test it is 2 wk. For the treatment described in Fig. 6, not only the immune system but also the tumor load returns to its original high level as illustrated by the treatment trajectory in the cancer–immunity phase space (Fig. 6*A*). The trajectory exemplifies the fact that although the radiation therapy reduces the tumor load (cancer level C) by 40%, after the treatment ends, C gradually increases (within a few weeks) back to its original (before the treatment) high level. Even if the radiation therapy is applied for 100 d instead of 30 d, so that C is further reduced by more than 50%, tumor load gradually returns to its original level (Fig. 7*A*, *Left*). It should be noted that a 100-d treatment is not realistic because it is usually above the maximum permissible dose for the radiation therapy.

### Translational Relevance.

The above results suggest that radiation therapy alone may not be effective. As we illustrate next, a likely more effective strategy would be based on alternating combined radiation and immune therapy in a two-stage approach—an effective therapy that can induce a transition from the high cancer state (H) to the intermediate state (I), H2I transitions, followed by an effective alternate therapy that can induce a transition from the intermediate cancer state (I) to the low cancer state (L), I2L transitions. Only in rare situations can direct H2L transitions be induced.

### Assessment of Simulated Immunotherapy.

The effect of hypothetical DC immunotherapy is shown in Fig. 7*A*, *Right*. In this hypothetical therapy, the C level decreases to only ∼85% of its original value at the end of therapy. However, because the therapy leads to a significant increase in the D level, it induces a transition from the (H) state to the (I) state (H2I transition; Fig. 7*A*, *Right*), though the induced transition is very slow. Inspection of the time dependence of the H2I trajectory in the phase plane reveals that it takes ∼3 mo for the system to approximate the unstable steady state (or the dynamical barrier) between the (H) and (I) states (hollow green circle in the phase plan). Upon crossing this dynamical barrier, the system gradually approaches the intermediate state and converges to this state after ∼2 y. Although the outcome of standalone immunotherapy is far better than that of the standalone radiation therapy, the slow convergence makes the process more sensitive to noise and increases the chance that random mutations or other events will cause the trajectory to change its course toward the (H) state. In addition, applications of too intense immunotherapy and/or application of immune boost for long periods of time increase the risk for side effects such as autoimmune diseases, making standalone strong immunotherapy less desirable.

### The Merit of Alternating Combined Treatments.

We proceed to test the efficiency of therapies that are based on alternating radiation therapy and immunotherapy. In the hypothetical example shown in Fig. 7*B*, *Left*, we simulate three consecutive alternating treatments. Each treatment consists of combining a set of 8 d of radiation therapy, followed by a set of 8 d of immunotherapy. As illustrated in Fig. 7*B*, *Left*, this hypothetical alternating combined therapy successfully induces H2I transition. Moreover, the alternating combined therapy uses far lower doses of radiation, and thus has significantly reduced side effects in comparison with a single standalone radiation therapy (60) and induces H2I transitions without the need to overboost the immune system (in comparison with a standalone immune therapy), which can have its own side effects.

Interestingly, inclusion of additional radiation therapy has an alarming effect—it causes a change in the therapy trajectory to converge back to the original (H) state (Fig. 7*B*, *Right*). Moreover, even if the total dose of radiation therapy is reduced, it still can cause the same negative effect (*SI Appendix*, Fig. S12*C*). We also test treatments in which both radiation therapy and immunotherapy are applied simultaneously on the same day, instead of consecutively on different days. In this case, the H2I transitions can be retained by using slightly reduced doses for radiation therapy and slightly milder amplification of the immune response (*SI Appendix*, Fig. S12 *A* and *B*). However, the potential benefits might not justify the practical difficulties of applying both therapies simultaneously.

### Therapeutic Protocols for Inducing I2L Transitions.

The dynamics of the cancer–immunity interplay in the vicinity of the (I) and (L) cancer states are very different from those in the vicinity of the (H) cancer state (Fig. 2*B*). For this reason, therapeutic protocols for inducing I2L transitions are expected to be different from those that are efficient for inducing H2I transitions. To devise an efficient protocol, we follow the same approach as before to assess the effectiveness of standalone radiation therapy, of standalone immunotherapy, and of alternating combined therapy. As illustrated in *SI Appendix*, Fig. S12*D*, standalone radiation therapy fails to induce the I2L transition similar to the inefficiency of radiation therapy to induce H2I transitions. However, standalone immunotherapy (five sets of 10 d of immunotherapy, with a 10-d break between any two consecutive sets) is found to be efficient to induce I2L transitions (Fig. 7*C*, *Left*). Surprisingly, alternating combined therapy turns out to be inefficient to induce the I2L (Fig. 7*C*, *Right*). The reason for the inefficiency is that the trajectory of the alternating combination is mostly orthogonal to the shortest path line that connects the (I) state with the dynamical barrier between the (I) and the (L) states (the unstable steady state marked by hollow green circle in Fig. 7*C*, *Right*).

### Effect of Time Delay on Treatment Efficacy.

To assess the effect of time delay on treatment efficacy, we convert the deterministic equations to include the time delay (*Methods*). Specifically, we assess the case of 15-d delay for exosome exchange between C and D. Treatment simulations (*SI Appendix*, Fig. S13) show that treatment trajectories are similar to those for the cases without considering time delay. However, most treatments become less effective, because time delay makes transitions harder (61), as discussed in the proceeding time delay section.

### Effects of Noise on Treatment Efficacy.

We also evaluate the effect of noise by turning the deterministic model into a stochastic model whose equations contain noise terms that represent white Gaussian noise (*Methods*). The presence of large noise can induce spontaneous transitions between the different quasi-stable cancer states. To test the noise effect on the dynamics itself, we select the noise level to be sufficiently low to avoid noise-induced spontaneous transitions during the duration of treatments. Still, noise affects treatment effectiveness. Starting with the example presented previously, in which a 25-d immunotherapy in the absence of noise can trigger transition from the (H) to the (I) state (Fig. 8*A*, solid green line), we find that the success rate for H2I transitions drops to 60% when the simulations start exactly from the levels of C, D, and K of the (H) state; the reason is that it is possible for the noise to divert the treatment trajectory away from the path to cross the dynamical barrier between the (H) and the (I) states (Fig. 8*A*, condition 3). However, in the presence of noise, the values of C, D, and K vary in the vicinity of the exact steady-state values before the treatment. Therefore, the initial point at which the therapy starts can be close to but not exactly at the (H) state. Hence, we continue to test the effects of the initial state on treatment outcome; to do so, we choose five different initial conditions 1–5 (Fig. 8*A*, *Insets*), whereby condition 3 is the exact (H) state, condition 1 is farthest from the transition state (the unstable steady state shown in hollow green circle) between the (H) and the (I) state, and condition 5 is closest to the dynamical barrier between the (H) and (I) states. As expected, the simulations reveal that the closer the initial state (when the therapy starts) to the dynamical barrier, the higher the probability to induce H2I transition, with the lowest probability being 20% when the therapy starts from condition 1. Notably, in the absence of noise, the therapy cannot induce H2I transition at all starting from this state. Considered together, the results indicate that noise sometimes can facilitate transitions during the treatment. In addition, when noise is included, there is low probability to induce direct H2L transitions that would not be possible in the absence of noise. Moreover, the results also demonstrate that it is crucial to start the treatment at the right time when the cancer–immunity interplay is at a state as close as possible to the dynamical barrier for transitions to be induced. It follows that the cancer and immune states should be monitored closely and that therapy should be individualized and timed to increase the likelihood of the treatment success.

## Conclusion and Discussion

In this study, we devised a theoretical framework for exosome exchange-based cancer–immunity interplay (the ECI model) to examine the contribution of exosome exchange to the intricate cancer–immunity interplay. We combined the complex interplay into a simplified generic model whose three components are cancer cells C, dendritic cells D (representing precursor, immature and mature types) and killer cells K (representing cytotoxic T cells, helper T cells, effector B cells, and natural killer cells). Distinct from other models (36, 38), the time evolution dynamics for each cell population was modeled with a set of nonlinear rate equations in a way similar to the analysis for gene regulatory circuits. Due to the specific cell–cell interactions among C, D, and K cells, the ECI model resembles the self-activating toggle switch, a gene circuit module that plays an important role in decision-making for many biological processes, e.g., epithelial–mesenchymal transitions (44, 62) and cell differentiation (63⇓–65). It has been shown that such a circuit topology allows for three stable steady states (44). For the current ECI model, we also found that the cancer–immune system, at a certain level of immune recognition *ρ* and when exosome exchange is included, can allow three cancer states: high (H), intermediate (I), and low (L). The model prediction regarding the (I) state is consistent with the immunoediting theory (45). Interestingly, owing to the special interactions among different cancer and immune cells, the immune system is most active in the intermediate state. Note that the ECI model does not particularly depend on the exosome exchange mechanism. Therefore, it also can be applied to many other mechanisms of cell–cell communications between the various cell types that are involved in the cancer–immunity interplay.

To study tumorigenesis, we incorporated the immune recognition represented by a variable *ρ*. During tumorigenesis, *ρ* gradually changes depending on the cancer stages. We proposed that, at the onset of a solid tumor, *ρ* starts from zero and increases to full recognition (*ρ* = 1) at a certain response rate. We found that, depending on how fast the DCs recognize the new tumor, cancer development ends up at high, intermediate, or low cancer states. Once the cancer population reaches a quasi-steady state, *ρ* gradually decreases due to immune escape. This dynamics of immune recognition are consistent with change in tumor immunogenicity in the immunoediting theory (45); however, it is still unclear if the equilibrium state exists, what derives and maintains the equilibrium, and how to detect the equilibrium state (5).

We further analyzed the multistability of the model with respect to *ρ* within 2D and 3D bifurcation diagrams, from which we identified four different phases. In each of the phases, the model can be at one state, or stay in one of the multiple states. It is worth noting that the current ECI model does not consider the situation in which the tumor actively escapes the immune suppression, which also contributes to reducing immune recognition. For example, the tumor can limit antigen presentation (66) during cancer treatment. We take into account these effects by modeling the time evolution of *ρ* in a future study.

We also performed stochastic analysis on the ECI model, from which we constructed a corresponding effective landscape. The effective landscape not only shows the relative stability of each steady state, but also provides an insight into cancer prognosis. Our findings imply that an effective treatment should not aim to simply reduce the size of the tumor, but to lead the dynamics of the cancer–immunity interplay to cross the dynamical barrier (effective potential barrier) during treatment from the (H) to the (I) state and from the (I) to the (L) states. In many cases, effective therapeutic protocols would be different for the H-to-I and the I-to-L transitions. Finally, we assessed the effects of time delay and biological noise on the design and effectiveness of therapeutic strategies. We note that noise effects of the cancer–immunity are still poorly understood (36).

The ECI model enables us to evaluate the effectiveness of various treatment strategies, based on the mathematical modeling of both the cancer–immunity interplay and the modeling of the response of cancer–immunity to the hypothetical therapy. Our current model suggests that standalone radiation therapy is mostly ineffective, as the system usually goes back to the original quasi-stable state after the treatment. Radiation therapy may reduce not only the tumor load but also the immune response (58, 59, 67), making the system lean more toward the basin of the high cancer state. We also found that the immunotherapy is more effective. However, we may overestimate the effectiveness of the immunotherapy due to the following two factors. First, the slow convergence of cancer–immunity interplay to the (L) state during a standalone immunotherapy may increase the risk of failure due to some unpredicted factors such as de novo mutations in the tumor. Second, we do not know the maximum tolerated dosage for the immunotherapy, so it is possible that too much drug was applied, causing intolerable side effects. For example, excessive use of immunotherapy may result in strong immune tolerance, and it can in turn reduce the immune recognition *ρ* or induce autoimmune diseases. Such potential effects will be taken into account in future studies.

Though the ECI model presented here turned out to be a valuable starting point, it calls for future studies in which the cell subpopulations in C, D, and K are individually included in the ECI model. For example, macrophages usually play a role in eliminating tumor cells and in priming naïve T and B cells by antigen presentation, whereas tumor-associated macrophages may stimulate tumor growth (10, 68). Moreover, regulatory T cells can induce immune tolerance, and may cause oscillatory cycles of the immune system (69). Presumably, the effectiveness of cancer treatment could be affected by the dynamics of the regulatory T cells; it is crucial to investigate their roles in the population dynamics of the cancer–immune system.

The ECI model has the potential to be valuable for clinical practice; for this to happen, one first has to adjust the model parameters for each individual patient. Only then, and in accordance with recommended treatments for the specific patient, can the model be used to assess possible optional combined strategies and propose the optimal one according to the patient-specific therapeutic trajectories as calculated by the ECI model. We hope that the predictive power of this approach of model-assisted personalized practice will be tested in the near future by dedicated trials.

## Methods

### Stability Analysis.

The dynamics of the ECI model is described by a set of 3D nonlinear rate equations in Eq. **1**. To compute the steady-state solutions, we calculate two nullclines—the first one satisfies dD/dt = 0 and dK/dt = 0 (solid navy line in Fig. 2), and the second one satisfies dC/dt = 0 and dK/dt = 0 (solid brown line in Fig. 2). The nullclines are projected onto the phase plane, constructed by the concentrations of the D (*x* axis) and C (*y* axis). The intersections of these two projected lines are the steady states for the whole system. The stability of the steady states can be further determined by the linear approximation method. The nullclines are computed by contour-based method, and the bifurcation diagrams are calculated by PITCON7 package (70).

### Simulation with Time Delay.

We consider the time delays caused by the exosome-mediated communications between the C and D cell populations. Instead of directly using Eq. **1**, we replace the terms

The delay differential equations are integrated by dde23 from MATLAB (71).

### Stochastic Simulation.

Eq. **1** can be rewritten in a vector form as

### Tumorigenesis and Treatment Modeling.

We specifically model the time evolution of the level of immune recognition *ρ* for the tumor onset process. *ρ* follows

where *ρ* initially starts at 0 and increases until it saturates at 1. The tumor onset is simulated by both Eq. **1** and Eq. **4**, and the initial condition for (C, D, K) is (0, 100, 0) (cells per microliter).

In the simulation of the radiation therapy, additional terms **1**, respectively, where **1**. The parameters are listed in *SI Appendix*, section 2.

## Acknowledgments

This work was supported by the Center for Theoretical Biological Physics sponsored by National Science Foundation Grants PHY-1427654 and NSF-MCB-1214457, the Cancer Prevention and Research Institute of Texas (CPRIT), the Tauber Family Funds, and the Magut-Glass Chair in Physics of Complex Systems at Tel-Aviv University. J.N.O. is a CPRIT Scholar in Cancer Research sponsored by CPRIT. M.L. is supported by a training fellowship from Keck Center for Interdisciplinary Bioscience Training of the Gulf Coast Consortia CPRIT Grant RP140113.

## Footnotes

↵

^{1}M.L. and B.H. contributed equally to this work.- ↵
^{2}To whom correspondence may be addressed. Email: jonuchic{at}rice.edu or eshelbj{at}gmail.com.

Author contributions: M.L., B.H., S.M.H., J.N.O., and E.B.-J. designed research; M.L. and B.H. performed research; M.L., B.H., S.M.H., J.N.O., and E.B.-J. analyzed data; and M.L., B.H., S.M.H., J.N.O., and E.B.-J. wrote the paper.

Reviewers: K.J.P., Johns Hopkins University; W.L., University of Maryland, College Park; and S.S., Hebrew University of Jerusalem.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1416745111/-/DCSupplemental.

## References

- ↵.
- Schlom J

- ↵
- ↵.
- Vande Woude GF,
- Klein G

- Mougiakakos D,
- Choudhury A,
- Lladser A,
- Kiessling R,
- Johansson CC

- ↵
- ↵
- ↵
- ↵
- ↵.
- Igney FH,
- Krammer PH

- ↵
- ↵
- ↵.
- Smyth MJ, et al.

- ↵.
- Nagorsen D,
- Scheibenbogen C,
- Marincola FM,
- Letsch A,
- Keilholz U

- ↵.
- Toes REM,
- Ossendorp F,
- Offringa R,
- Melief CJM

- ↵.
- Kuby J

- ↵
- ↵
- ↵.
- Schroder K,
- Hertzog PJ,
- Ravasi T,
- Hume DA

- ↵
- ↵
- ↵
- ↵
- ↵.
- Vaupel P,
- Kallinowski F,
- Okunieff P

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Zhang H-G,
- Grizzle WE

- ↵
- ↵
- ↵
- ↵
- ↵.
- Yu S, et al.

- ↵.
- Liu C, et al.

- ↵.
- Capasso V,
- Lachowicz M

- Banasiak J,
- Chaplain MAJ,
- Miękisz J

- ↵
- ↵
- ↵
- ↵.
- Pappalardo F,
- Pennisi M,
- Ricupito A,
- Topputo F,
- Bellone M

- ↵.
- Wilkie KP,
- Hahnfeldt P

- ↵
- ↵
- ↵
- ↵.
- Lu M, et al.

- ↵
- ↵
- ↵.
- Gerosa F, et al.

- ↵.
- Vivier E, et al.

- ↵.
- Nelson BH

- ↵.
- Clayton A,
- Mitchell JP,
- Court J,
- Mason MD,
- Tabi Z

- ↵
- ↵
- ↵
- ↵.
- Menetrier-Caux C, et al.

- ↵.
- Scarlett UK, et al.

- ↵.
- Wang J,
- Zhang K,
- Xu L,
- Wang E

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Lu M,
- Jolly MK,
- Levine H,
- Onuchic JN,
- Ben-Jacob E

- ↵
- ↵.
- Zhang B,
- Wolynes PG

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Kloeden PE,
- Platen E

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Biophysics and Computational Biology