Towards a unified theory of cancer risk

Martin Nowak and Bartlomiej Waclaw conclude their recent commentary [1] on the “bad luck and cancer” debate with a look to the future:

“The earlier analysis by Tomasetti and Vogelstein has already stimulated much discussion… It will take many years to answer in detail the interesting and exciting questions that have been raised.”

I agree. When a couple of journalists [2, 3] contacted me for comments on the latest follow-up paper from Christian Tomasetti, Bert Vogelstein and Lu Li, I emphasized what can be gained from rekindling the decades-old debate about the contribution of extrinsic (or environmental, or preventable) factors to cancer risk. In particular, the diverse scientific critiques of Tomasetti and Vogelstein’s analysis suggest important avenues for further inquiry.

My own take is summarized in the figure below. This diagram (inspired by Tinnbergen) reframes the question in terms of proximate mechanisms and ultimate causes. It also provides a way of categorizing cancer etiology research.

Causes of cancer

Tomasetti and Vogelstein’s 2015 paper [4] demonstrated that the lifetime number of stem cell divisions is correlated with cancer risk across human tissues (part A in the figure). Colleagues and I have argued [5, 6] that, although characterizing this association is important, it cannot be used to infer what proportion of cancer risk is due to intrinsic versus extrinsic factors. This is because cancer initiation depends not only on mutated cells, but also on the fitness landscape that governs their fate, which is determined by a microenvironment that differs between tissues (figure part B).

Moreover, the supply of mutated cells and the microenvironment are both shaped by an interaction of nature and nurture (figure part C). In a recently published paper [7], Michael Hochberg and I draw attention to the relationship between cancer incidence and environmental changes that alter organism body size and/or life span, disrupt processes within the organism, or affect the germline (figure part D). We posit that “most modern-day cancer in animals – and humans in particular – are due to environments deviating from central tendencies of distributions that have prevailed during cancer resistance evolution”. We support this claim in our paper with a literature survey of cancer across the tree of life, and with an estimate of cancer incidence in ancient humans based on mathematical modelling [7].

To understand why cancer persists at a certain baseline level even in stable environments, we must further examine the role of organismal evolution (figure part E). If cancer lowers organismal fitness then we might expect selection for traits that reduce risk. But continual improvement in cancer prevention is expected to come at a cost, and the net effect on fitness will depend on life history. For example, more stringent control of cell proliferation might reduce cancer risk and so lower the mortality rate at older ages, while also increasing deaths in juveniles and young adults due to impaired wound healing. We can predict outcomes of such trade-offs by calculating selection gradients, which is what I’ve been doing in a research project that I presented at an especially stimulating MBE conference in the UK last week.

The quest to understand cancer risk must then encompass not only cell biology, but also ecology and evolution at both tissue and organismal levels. One of my goals is to make connections between these currently disparate lines of research in pursuit of a more unified theory.


  1. Nowak, M. A., & Waclaw, B. (2017). Genes, environment, and “bad luck”. Science, 355(6331), 1266–1267.
  2. Ledford, H. (2017) DNA typos to blame for most cancer mutationsNature News.
  3. Chivers, T. (2017) Here’s Why The “Cancer Is Caused By Bad Luck” Study Isn’t All It Seems. Buzzfeed.
  4. Tomasetti, C., & Vogelstein, B. (2015). Variation in cancer risk among tissues can be explained by the number of stem cell divisions. Science, 346(6217), 78–81.
  5. Noble, R., Kaltz, O., & Hochberg, M. E. (2015). Peto’s paradox and human cancers. Philosophical Transactions of the Royal Society B: Biological Sciences, 370(1673), 20150104–20150104.
  6. Noble, R., Kaltz, O., Nunney, L., & Hochberg, M. E. (2016). Overestimating the Role of Environment in Cancers. Cancer Prevention Research, 9(10), 773–776.
  7. Hochberg, M. E., & Noble, R. J. (2017). A framework for how environment contributes to cancer risk. Ecology Letters20(2), 117–134.

The Box-Einstein surface of mathematical models

As a mathematical modeller in evolutionary biology, my seminar bingo card has four prime boxes. Watching a talk about evolution, I count down the minutes to the first appearance of Dobzhansky’s “nothing in biology” quote (or some variant thereof) or a picture of Darwin’s “I think” sketch. For mathematical modelling, it’ll be either Albert Einstein or George Box:

“All models are wrong but some are useful” – George Box

“Everything should be made as simple as possible, but not simpler” – probably not Albert Einstein

Of course, such quotes are popular for good reason, and I’m not criticising those who use them to good effect, but all the same it can be fun to try to find a new way of presenting familiar material. That’s why in spring 2015 I came up with and tweeted a visual summary of the latter two aphorisms, which I named the Box-Einstein surface of mathematical models:


The grey region in the plot ensures that all possible models have some degree of “wrongness”, but the contours in the remaining region tell us that some models are useful all the same. To find the most useful description of a particular phenomenon, we must reduce complexity without overly increasing wrongness.

A key thing to understand about this diagram is that although the boundary of the grey region is invariant, the surface is changeable. If our empirical knowledge of the system becomes richer, or if we change the scope of our enquiry, the most useful model may be more or less complex than before.

Einstein’s quote can be seen as simply paraphrasing Occam’s razor, but I think it has additional meaning with regard to what Artem Kaznatcheev calls heuristic and abstract mathematical models, such as are generally used in biology. In statistics, a simple model has few degrees of freedom, which is desirable to reduce overfitting. However, statisticians should also beware what JP Simmons and colleagues termed “researcher degrees of freedom”:

“In the course of collecting and analyzing data, researchers have many decisions to make: Should more data be collected? Should some observations be excluded? Which conditions should be combined and which ones compared? Which control variables should be considered? Should specific measures be combined or transformed or both?

“It is rare, and sometimes impractical, for researchers to make all these decisions beforehand. Rather, it is common (and accepted practice) for researchers to explore various analytic alternatives, to search for a combination that yields “statistical significance,” and to then report only what “worked.” The problem, of course, is that the likelihood of at least one (of many) analyses producing a falsely positive finding at the 5% level is necessarily greater than 5%.”

Likewise, when a researcher makes a mathematical model of a dynamical system – be it a set of differential equations or a stochastic agent-based model – he or she makes numerous decisions, usually with more or less full knowledge of the empirical data against which the model will be judged.

But there’s an important difference between the process of collecting data and that of creating a mathematical model. Ideally, the experimentalist can minimise researcher degrees of freedom by following a suitable experimental design and running controls that enable him or her to test a hypothesis against a null according to a predetermined statistical model. For most mathematical models there is no such template, and a process of trial and improvement is unavoidable, forgivable, and even desirable (inasmuch as it strengthens understanding of why the model works). The role of mathematical modeller is somewhere between experimentalist and pure mathematician. By making our models as simple as possible, we shift ourselves further toward the latter role, and our experimentation becomes less about exploiting our freedom and more about honing our argument.

For further reading, check out Artem Kaznatcheev’s insightful post about what “wrong” might mean, and why Box’s quote doesn’t necessarily apply to all types of model.

Visualizing evolutionary dynamics with ggmuller

A how-to guide for making beautiful plots of your data.

Readers of experimental evolution studies by Richard Lenski, Jeffrey Barrick and others, or of cancer evolution reviews, will have seen stacked area plots that combine information about frequency dynamics and phylogeny. The example below is from a recent preprint by Rohan Maddamsetti, Lenski and Barrick. The horizontal axis is time (measured in generations), and each coloured shape represents a genotype, with height corresponding to genotype frequency (or relative abundance) at the corresponding time point. Descendant genotypes are shown emerging from inside their parents.


Such diagrams are sometimes termed Muller plots in honour of Hermann Joseph Muller (of ratchet fame), who used them in 1932 to illustrate an evolutionary advantage of sex. I wanted to draw Muller plots of some evolutionary simulations. Unable to find generic software, I spent a few days coding in R and tweeted an example result:


Encouraged by the response, I expanded my code into the fully documented R package ggmuller (thanks to Sean Leonard for suggesting the name). Ggmuller is my attempt to set an easy-to-use standard for representing evolutionary dynamics. You can install the package from github with a few lines of R code:

install.packages("devtools") # enables installation from github

Some of the code that follows depends on the latest version of ggmuller; if you installed the package before 20th August then please reinstall it for full functionality.

How does it work?

The core function is get_Muller_df, which constructs a specially ordered data frame from which we can draw a stacked area plot. get_Muller_df takes two data frames as input: first an adjacency matrix that defines a phylogeny (with columns named “Identity” and “Parent”), and second the populations of each genotype over time (with columns named “Identity”, “Generation” and “Population”). The genotype names in the “Identity” and “Parent” columns can be any words or numbers. You can use a phylo object instead of an adjacency matrix, provided the genotype names in the population data are integers that correspond to the node numbers in the phylo object.

get_Muller_df records a path through the phylogenetic tree. Each genotype appears exactly twice in the path: once before and once after all of the genotype’s descendants. The function then associates each instance of each genotype with half of the genotype’s frequency over time.

To draw plots, ggmuller exploits Hadley Wickham’s ggplot package (which accounts for the first part of its name). Because ggplot’s stacked area plot respects row order, and because the two instances of each genotypes are coloured identically, the resulting plot shows each genotype emerging from the centre of its parent (following Jeffrey Barrick’s example, multiple descendants are stacked from top to bottom in chronological order of appearance). Whereas I could have made descendants emerge at different heights – as in the Lenski lab figure above – it seems simpler and less ambiguous to have all descendants emerge exactly in the middle of their parents. I propose this as a standard to facilitate comparing plots of this type.

Basic usage

If your adjacency matrix and your population data frame are properly formatted, you can use ggmuller to visualize your results with just two lines of R code:

Muller_df <- get_Muller_df(example_edges, example_pop_df)

Simply replace example_edges and example_pop_df with the names of the data frames that contain your adjacency matrix and population data, respectively.

That’s all you need for basic usage. The rest of this post explains why these plots are such a powerful way of visualizing evolutionary dynamics, and introduces some additional features and options.

A couple of minimal examples

To get a better idea of why you might want to use ggmuller, let’s recreate Figure 1 of a recent commentary by Mark Robertson-Tessi & Sandy Anderson. The figure (drawn by Chandler Gatenbee) comprises two diagrams contrasting the “traditional” selective sweep model of cancer progression and the“Big Bang” neutral evolution model proposed by Andrea Sottoriva and colleagues.

Most of the code that follows is to create an appropriate data set. First we make an adjacency matrix that defines the phylogeny. In this minimal example, the gentoype names are integers, and each genotype has exactly one descendant (1 begets 2, which begets 3, which begets 4):

edges1 <- data.frame(Parent = 1:3, Identity = 2:4)

Next we construct a data frame of genotype populations over time. In the neutral evolution case, all genotypes (except the background genotype) have the same growth rate:

# a function for generating exponential growth curves:
pop_seq <- function(gens, lambda, start_gen) c(rep(0, start_gen),
 exp(lambda * gens[0:(length(gens) - start_gen)]))
gens <- 0:150 # generations
lambda <- 0.1 # baseline exponential growth rate
pop1 <- data.frame(Generation = rep(gens, 4),
 Identity = rep(1:4, each = length(gens)), 
 Population = c(1E2 * pop_seq(gens, lambda, 0), 
 pop_seq(gens, 2*lambda, 0), 
 pop_seq(gens, 2*lambda, 3), 
 pop_seq(gens, 2*lambda, 8)))

We use get_Muller_df to combine the phylogeny and population dynamics into a data frame we can plot from:

Muller_df1 <- get_Muller_df(edges1, pop1)

We create the plot using Muller_plot (which is a wrapper for ggplot’s stacked area plot) using a custom colour palette (optional):

my_palette <- c("white", "skyblue2", "red", "green3")
plot1 <- Muller_plot(Muller_df1, palette = my_palette)

We proceed similarly for the second panel, in which genotypes have different growth rates (for example, due to the accumulation of advantageous mutations):

edges2 <- data.frame(Parent = 1:3, 
 Identity = 2:4)
pop2 <- data.frame(Generation = rep(gens, 4),
 Identity = rep(1:4, each = length(gens)),
 Population = c(1E2 * pop_seq(gens, lambda, 0), 
 pop_seq(gens, 2*lambda, 0), 
 pop_seq(gens, 4*lambda, 40), 
 pop_seq(gens, 8*lambda, 80)))
Muller_df2 <- get_Muller_df(edges2, pop2)
plot2 <- Muller_plot(Muller_df2, palette = my_palette)

Finally, we draw the plots together using the gridExtra package:

grid.arrange(plot1, plot2)


And here we see the power of the Muller vizualisation: even though these systems have identical phylogenies, we can see at a glance that they have very different evolutionary dynamics.

A more sophisticated phylogeny

Now let’s look at a more interesting, branched phylogeny:

edges3 <- data.frame(Parent = paste0("clone_", 
 LETTERS[c(rep(1:3, each = 2), 2, 5)]), 
 Identity = paste0("clone_", LETTERS[2:9]))

The population data are generated using diverse fitness values:

fitnesses <- c(1, 2, 2.2, 2.5, 3, 3.2, 3.5, 3.5, 3.8)
pop3 <- data.frame(Generation = rep(gens, 9),
 Identity = paste0("clone_", LETTERS[rep(1:9, each = length(gens))]),
 Population = c(1E2 * pop_seq(gens, fitnesses[1]*lambda, 0), 
 pop_seq(gens, fitnesses[2]*lambda, 0), 
 pop_seq(gens, fitnesses[3]*lambda, 10), 
 pop_seq(gens, fitnesses[4]*lambda, 20),
 pop_seq(gens, fitnesses[5]*lambda, 30),
 pop_seq(gens, fitnesses[6]*lambda, 40),
 pop_seq(gens, fitnesses[7]*lambda, 50),
 pop_seq(gens, fitnesses[8]*lambda, 50),
 pop_seq(gens, fitnesses[9]*lambda, 60)),
 Fitness = rep(fitnesses, each = length(gens)))

Of course there are more efficient ways of coding this in R, but perhaps the long-winded version is clearer for readers unfamiliar with lapply and its siblings. Note the inclusion of an optional “Fitness” column containing the fitness values. We combine the data using get_Muller_df as before:

Muller_df3 <- get_Muller_df(edges3, pop3)

Finally, we create two versions of the plot, first coloured by genotype identity (the default) and then (using the RColorBrewer package) coloured by fitness, making use of that optional column:

plot3 <- Muller_plot(Muller_df3, add_legend = TRUE)
num_cols <- length(unique(Muller_df3$Fitness)) + 1
Muller_df3$Fitness <- as.factor(Muller_df3$Fitness)
plot3a <- Muller_plot(Muller_df3, colour_by = "Fitness", 
 palette = rev(colorRampPalette(brewer.pal(9, "YlOrRd"))(num_cols)), 
 add_legend = TRUE)
grid.arrange(plot3, plot3a)


 Optional extras

Since simulations can result in very large sets of genotypes, many of which never become abundant, get_Muller_df has a threshold option to exclude rarities. By default, the data frame includes only rows with nonzero population; this can be overridden by setting add_zeroes = TRUE, in which case all genotypes are represented at every generation (which makes it easier to add new columns, for example).

Converting between adjacency matrix and phylo representations

We can visualize the tree from our branched phylogeny example using Emmanuel Paradis’ ape package, but first we need to use the ggmuller function adj_matrix_to_tree to convert our adjacency matrix (with arbitrary node names) to a phylo object that ape can understand. Phylo is the standard way of representing trees in R, and it requires nodes to be numbered in a particular order.

tree <- adj_matrix_to_tree(edges3)

After optionally labeling the nodes, we use ape’s plot.phylo function to draw the tree:

tree$tip.label <- 1:length(tree$tip.label) # optional
tree$node.label <- (length(tree$tip.label) + 1):10 # optional
plot(tree, show.node.label = TRUE, show.tip.label = TRUE, tip.color = "red")


If you already have a tree in phylo format, you can use that instead of an adjacency matrix as input to get_Muller_df, as long as the genotype names in the population data are integers that correspond to the node numbers in the phylo object.

Showing changes in population size

Although ggmuller deals in frequencies, you can hack it to show changes in the total population size by adding a dummy common ancestor to your phylogeny, with population equal to the maximum total population size minus the total population size at each time point, and then putting a new scale on the y-axis. By default, the dummy common ancestor (which corresponds to unfilled space) will be shown in black, whereas actual genotypes will have paler colours. I may automate this process in future versions.

Next steps

I hope ggmuller will prove useful to researchers in experimental evolution, cancer research, population genetics, and beyond. I plan to submit the package to CRAN and will continue to add functionality. Meanwhile, if you use ggmuller in a publication, please cite the github version including doi:10.5281/zenodo.240589, and please notify me via email or twitter so I can keep track of uptake. Bug reports and feature requests would be most welcome.

[Edited on 5th September 2016 to fix a typo in the code defining “edges3”.]