Dark selection from spatial cytokine signaling networks

Greetings, Theory, Evolution, and Games Group! It’s a pleasure to be on the other side of the keyboard today. Many thanks to Artem for the invite to write about some of our recent work and the opportunity to introduce myself via this post. I do a bit of blogging of my own over at vcannataro.com — mostly about neat science I stumble over while figuring out my way.

I’m a biologist. I study the evolutionary dynamics within somatic tissue, or, how mutations occur, compete, accumulate, and persist in our tissues, and how these dynamics manifest as aging and cancer (Cannataro et al., 2017a). I also study the evolutionary dynamics within tumors, and the evolution of resistance to targeted therapy (Cannataro et al., 2017b).

In November 2016 I attended the Integrated Mathematical Oncology Workshop on resistance, a workweek-long intensive competitive workshop where winners receive hard-earned $$ for research, and found myself placed in #teamOrange along with Artem. In my experience at said workshop (attended 2015 and 2016), things usually pan out like this: teams of a dozen or so members are assembled by the workshop organizers, insuring a healthy mix of background-education heterogeneity among groups, and then after the groups decide on a project they devise distinct but intersecting approaches to tackle the problem at hand. I bounced around a bit early on within #teamOrange contributing to our project where I could, and when the need for a spatially explicit model of cytokine diffusion and cell response came up I jumped at the opportunity to lead that endeavor. I had created spatially explicit cellular models before — such as a model of cell replacement in the intestinal crypt (Cannataro et al., 2016) — but never one that incorporated the diffusion or spread of some agent through the space. That seemed like a pretty nifty tool to add to my research kit. Fortunately, computational modeler extraordinaire David Basanta was on our team to teach me about modeling diffusion (thanks David!).

Below is a short overview of the model we devised.

Let me briefly motivate this work. For more details, see our preprint (Kaznatcheev et al., 2017) — especially section 3.3 — or Artem’s first post on this topic. Ruxolitinib therapy resistance in chronic myelomonocytic leukemia (CMML) does not follow the often-assumed trajectory of resistance—that of a subclone with a genetic variant conveying resistance outcompeting nonresistant clones in the presence of treatment. In our system, resistance arises, but with no measurable change in genetic profile or drastic change in population size. The selection for resistance is hidden, hence the “dark selection” theme of our project (and we really went to town on the “dark side of selection” analogies in our final presentation). So, how does resistance arise? We already discussed three potential mechanisms, and I want to focus on my favorite.

We know that the therapy disrupts cytokine signaling, and that signaling pathways and receptors (especially in regards to JAK/STAT) are dynamic over time even in the absence of genetic alteration, e.g. cells may become ‘desensitized’ to a molecule that is oversaturated or hypersensitive to a molecule that is largely absent. One plausible mechanism of resistance, we think, is that this therapy may be selecting for resistance at the phenotypic level. Perhaps the drug prevents cytokine signaling and cell activation via a default pathway, but the absence of signaling over time creates hypersensitive cells that may become ‘active’ through an alternative pathway (again, please see the preprint for biological details! Or Robert’s post on potential molecular implementations). Perhaps a phenotypic switch occurs.

As a preliminary step (and with one day left in the competition) we can use an in silico model to investigate the plausibility of this hypothesis^[1] and generate testable hypotheses about different mechanisms by which this phenotypic switch occurs. For instance, if the probability of the phenotypic switch was a function of the local cytokine load, how would the resultant dynamics manifest differently?

So, along with the other approaches, we modeled non-moving immortal cells arranged in space that release and respond to cytokines. Cells are ‘active’ and release cytokines, or if local cytokine levels become too high cells become ‘inactive’ and stop releasing cytokines for a period of time. Inactive cells, after a rest period, become activated once their local cytokines levels surpass some specified activation threshold. All the while, cytokines are diffusing through this space.^[2]

The resultant dynamics are really interesting. Check out the video we showed during our final presentation. Each box represents a 75×75 grid of cells, with the left box depicting the activation status of the cells (black=active, gray=inactive) and the right box depicting cytokine load (red=high, white=low). A random distribution of cells start active, and cytokines spread through the system. Cells become active and inactive depending on the cytokine levels, with the dynamics taking on a behavior reminiscent of Conway’s game of life. Eventually, an cyclic equilibrium is reached.

During the 26th second of the video, we introduce a drug that turns off the activated status of the cells and the ability of the cells to become activated by cytokines utilizing the previous cellular pathway. We see cytokine level decay. At the same time, the cells are able to phenotype switch to being activated by cytokines via a different pathway, and become active at some lower threshold of cytokines. This phenotypic switching occurs randomly throughout the simulation, and eventually cytokine levels return to pre-drug levels. You can see this below.

We model this rate of switching as either local cytokine level independent or occurring at a rate proportional to local cytokine levels. We see a more rapid return to pre-drug cytokine levels for the latter scenario. Furthermore, the cyclic pattern of cytokine equilibrium has a smaller range when compared to pre-drug levels. This is more drastic with the low spread in the cell-autonomous condition (left panel) versus the moderate amount of variance in final cytokine activity in the non-cell autonomous condition (right panel).^[3]

Overall, we showed that a symptom of CMML (high cytokine loads) may decline and then return of pre-drug levels in the presence of drug through a biological plausible mechanism of phenotypic switching, even in simulations with constant cell population sizes. I hope that this hypothesis, along with others interrogated with mathematical modeling during the IMO workshop, will be tested soon in vivo. The results of these experiments will expand our understanding of cancer growth, drug resistance, and evolutionary biology.

One last note: If you happen to find yourself in the position to attend an IMO workshop, I highly recommend applying for a travel award and attending. Especially as a graduate student or postdoc. Not only can you learn new skills (as I detailed above), you will also meet and collaborate with an international and diverse group of researchers. And, you may even get a preprint or publication out of the experience!

Notes and References

1. Given some time and $, I think that we can test the cytokine-network hypothesis in vivo (or at least in vitro), and that is our eventual plan.
2. The trick to adding diffusion was relatively straightforward. Our simulation takes place on a two dimensional grid, and every space in the grid has some level of cytokines. At the next iteration in time, for each space, cytokines enter that space from the surrounding eight spaces at some rate times the amount of cytokine in every surrounding space, and leave the space at eight times the rate of cytokine spread times the amount within the space. Thus, cytokines diffuse to all surrounding space while also entering from all surrounding space. On top of this, cytokines also naturally decay at some specified rate after spreading.
3. This second observation of change in range is reminiscent of the difference in variance that Artem predicts from a non-spatial model (red vs. blue lines in the final figure of the three mechanisms post). However, the first observation of faster response in the non-cell autonomous case is harder to reconcile, since Artem specifically selected activation functions for cell-autonomous vs non-cell-autonomous that would produce the same response time.

Cannataro, V. L., McKinley, S. A., & St Mary, C. M. (2016). The implications of small stem cell niche sizes and the distribution of fitness effects of new mutations in aging and tumorigenesis. Evolutionary Applications, 9(4), 565-582.

Cannataro, V. L., McKinley, S. A., & St Mary, C. M. (2017a). The evolutionary trade‐off between stem cell niche size, aging, and tumorigenesis. Evolutionary Applications.

Cannataro, V. L., Gaffney, S. G., Stender, C., Zhao, Z. M., Philips, M., Greenstein, A., & Townsend, J. P. (2017b). The likelihood of heterogeneity or additional mutation in KRAS or associated oncogenes to compromise targeting of oncogenic KRAS G12C. bioRxiv, 149724.

Kaznatcheev, A., Grimes, D. R., Vander Velde, R., Cannataro, V. L., Baratchart, E., Dhawan, A., … & Basanta, D. (2017). Dark selection for JAK/STAT-inhibitor resistance in chronic myelomonocytic leukemia. bioRxiv, 211151.