Mutation-selection networks of cancer initiation: tumor suppressor genes and chromosomal instability

Graph-structured populations elucidate the role of deleterious mutations in long-term evolution

Birth-death models are used to understand the interplay of genetic drift and natural selection. While well-mixed populations remain unaffected by the order of birth and death and where selection acts, evolutionary outcomes in spatially structured populations are affected by these choices. We show that the choice of individual moving to vacant sites-parent or offspring-controls the initial mutant placement on a graph and hence alters its fixation probability. Moving parent individuals introduces, to our knowledge, previously unexplored update rules and fixation categories for heterogeneous graphs. We identify a class of graphs, amplifiers of fixation, where fixation probability is larger than in well-mixed populations, regardless of the mutant fitness. Under death-Birth parent moving, the star graph is an amplifier of fixation, with a non-zero fixation probability for deleterious mutants, in contrast to very large well-mixed populations. Most Erdős-Rényi graphs of size 8 are amplifiers of fixation under death-Birth parent moving, but suppressors of fixation under Birth-death offspring moving. Surprisingly, amplifiers of fixation attain lower fitness in long-term evolution, despite favouring beneficial mutants, while suppressors of fixation attain higher fitness. These counterintuitive findings are explained by the fate of deleterious mutations and highlight the crucial role of deleterious mutants for adaptive evolution.

Evolutionary graph theory beyond single mutation dynamics: on how network-structured populations cross fitness landscapes

Spatially resolved datasets are revolutionizing knowledge in molecular biology, yet are under-utilized for questions in evolutionary biology. To gain insight from these large-scale datasets of spatial organization, we need mathematical representations and modeling techniques that can both capture their complexity, but also allow for mathematical tractability. Evolutionary graph theory utilizes the mathematical representation of networks as a proxy for heterogeneous population structure and has started to reshape our understanding of how spatial structure can direct evolutionary dynamics. However, previous results are derived for the case of a single new mutation appearing in the population and the role of network structure in shaping fitness landscape crossing is still poorly understood. Here we study how network-structured populations cross fitness landscapes and show that even a simple extension to a two-mutational landscape can exhibit complex evolutionary dynamics that cannot be predicted using previous single-mutation results. We show how our results can be intuitively understood through the lens of how the two main evolutionary properties of a network, the amplification and acceleration factors, change the expected fate of the intermediate mutant in the population and further discuss how to link these models to spatially resolved datasets of cellular organization.

Self-replicating artificial neural networks give rise to universal evolutionary dynamics

In evolutionary models, mutations are exogenously introduced by the modeler, rather than endogenously introduced by the replicator itself. We present a new deep-learning based computational model, the self-replicating artificial neural network (SeRANN). We train it to (i) copy its own genotype, like a biological organism, which introduces endogenous spontaneous mutations; and (ii) simultaneously perform a classification task that determines its fertility. Evolving 1,000 SeRANNs for 6,000 generations, we observed various evolutionary phenomena such as adaptation, clonal interference, epistasis, and evolution of both the mutation rate and the distribution of fitness effects of new mutations. Our results demonstrate that universal evolutionary phenomena can naturally emerge in a self-replicator model when both selection and mutation are implicit and endogenous. We therefore suggest that SeRANN can be applied to explore and test various evolutionary dynamics and hypotheses.

Aneuploidy Can Be an Evolutionary Diversion on the Path to Adaptation

Aneuploidy is common in eukaryotes, often leading to decreased fitness. However, evidence from fungi and human tumur cells suggests that specific aneuploidies can be beneficial under stressful conditions and facilitate adaptation. In a previous evolutionary experiment with yeast, populations evolving under heat stress became aneuploid, only to later revert to euploidy after beneficial mutations accumulated. It was therefore suggested that aneuploidy is a "stepping stone" on the path to adaptation. Here, we test this hypothesis. We use Bayesian inference to fit an evolutionary model with both aneuploidy and mutation to the experimental results. We then predict the genotype frequency dynamics during the experiment, demonstrating that most of the evolved euploid population likely did not descend from aneuploid cells, but rather from the euploid wild-type population. Our model shows how the beneficial mutation supply-the product of population size and beneficial mutation rate-determines the evolutionary dynamics: with low supply, much of the evolved population descends from aneuploid cells; but with high supply, beneficial mutations are generated fast enough to outcompete aneuploidy due to its inherent fitness cost. Our results suggest that despite its potential fitness benefits under stress, aneuploidy can be an evolutionary "diversion" rather than a "stepping stone": it can delay, rather than facilitate, the adaptation of the population, and cells that become aneuploid may leave less descendants compared to cells that remain diploid.

Mutation pressure, drift, and the pace of molecular coevolution

Most aspects of the molecular biology of cells involve tightly coordinated intermolecular interactions requiring specific recognition at the nucleotide and/or amino acid levels. This has led to long-standing interest in the degree to which constraints on interacting molecules result in conserved vs. accelerated rates of sequence evolution, with arguments commonly being made that molecular coevolution can proceed at rates exceeding the neutral expectation. Here, a fairly general model is introduced to evaluate the degree to which the rate of evolution at functionally interacting sites is influenced by effective population sizes (Ne), mutation rates, strength of selection, and the magnitude of recombination between sites. This theory is of particular relevance to matters associated with interactions between organelle- and nuclear-encoded proteins, as the two genomic environments often exhibit dramatic differences in the power of mutation and drift. Although genes within low Ne environments can drive the rate of evolution of partner genes experiencing higher Ne, rates exceeding the neutral expectation require that the former also have an elevated mutation rate. Testable predictions, some counterintuitive, are presented on how patterns of coevolutionary rates should depend on the relative intensities of drift, selection, and mutation.

Spectral dynamics of guided edge removals and identifying transient amplifiers for death-Birth updating

The paper deals with two interrelated topics: (1) identifying transient amplifiers in an iterative process, and (2) analyzing the process by its spectral dynamics, which is the change in the graph spectra by edge manipulation. Transient amplifiers are networks representing population structures which shift the balance between natural selection and random drift. Thus, amplifiers are highly relevant for understanding the relationships between spatial structures and evolutionary dynamics. We study an iterative procedure to identify transient amplifiers for death-Birth updating. The algorithm starts with a regular input graph and iteratively removes edges until desired structures are achieved. Thus, a sequence of candidate graphs is obtained. The edge removals are guided by quantities derived from the sequence of candidate graphs. Moreover, we are interested in the Laplacian spectra of the candidate graphs and analyze the iterative process by its spectral dynamics. The results show that although transient amplifiers for death-Birth updating are generally rare, a substantial number of them can be obtained by the proposed procedure. The graphs identified share structural properties and have some similarity to dumbbell and barbell graphs. We analyze amplification properties of these graphs and also two more families of bell-like graphs and show that further transient amplifiers for death-Birth updating can be found. Finally, it is demonstrated that the spectral dynamics possesses characteristic features useful for deducing links between structural and spectral properties. These feature can also be taken for distinguishing transient amplifiers among evolutionary graphs in general.

A spatial mutation model with increasing mutation rates

We consider a spatial model of cancer in which cells are points on the d-dimensional torus $\mathcal{T}=[0,L]^d$ , and each cell with $k-1$ mutations acquires a kth mutation at rate $\mu_k$ . We assume that the mutation rates $\mu_k$ are increasing, and we find the asymptotic waiting time for the first cell to acquire k mutations as the torus volume tends to infinity. This paper generalizes results on waiting for $k\geq 3$ mutations in Foo et al. (2020), which considered the case in which all of the mutation rates $\mu_k$ are the same. In addition, we find the limiting distribution of the spatial distances between mutations for certain values of the mutation rates.

Categorizing update mechanisms for graph-structured metapopulations

The structure of a population strongly influences its evolutionary dynamics. In various settings ranging from biology to social systems, individuals tend to interact more often with those present in their proximity and rarely with those far away. A common approach to model the structure of a population is evolutionary graph theory. In this framework, each graph node is occupied by a reproducing individual. The links connect these individuals to their neighbours. The offspring can be placed on neighbouring nodes, replacing the neighbours-or the progeny of its neighbours can replace a node during the course of ongoing evolutionary dynamics. Extending this theory by replacing single individuals with subpopulations at nodes yields a graph-structured metapopulation. The dynamics between the different local subpopulations is set by an update mechanism. There are many such update mechanisms. Here, we classify update mechanisms for structured metapopulations, which allows to find commonalities between past work and illustrate directions for further research and current gaps of investigation.

Simple mechanisms for the evolution of protein complexity

Proteins are tiny models of biological complexity: specific interactions among their many amino acids cause proteins to fold into elaborate structures, assemble with other proteins into higher-order complexes, and change their functions and structures upon binding other molecules. These complex features are classically thought to evolve via long and gradual trajectories driven by persistent natural selection. But a growing body of evidence from biochemistry, protein engineering, and molecular evolution shows that naturally occurring proteins often exist at or near the genetic edge of multimerization, allostery, and even new folds, so just one or a few mutations can trigger acquisition of these properties. These sudden transitions can occur because many of the physical properties that underlie these features are present in simpler proteins as fortuitous by-products of their architecture. Moreover, complex features of proteins can be encoded by huge arrays of sequences, so they are accessible from many different starting points via many possible paths. Because the bridges to these features are both short and numerous, random chance can join selection as a key factor in explaining the evolution of molecular complexity.

Optimal pathways control fixation of multiple mutations during cancer initiation

Cancer starts after initially healthy tissue cells accumulate several specific mutations or other genetic alterations. The dynamics of tumor formation is a very complex phenomenon due to multiple involved biochemical and biophysical processes. It leads to a very large number of possible pathways on the road to final fixation of all mutations that marks the beginning of the cancer, complicating the understanding of microscopic mechanisms of tumor formation. We present a new theoretical framework of analyzing the cancer initiation dynamics by exploring the properties of effective free-energy landscape of the process. It is argued that although there are many possible pathways for the fixation of all mutations in the system, there are only a few dominating pathways on the road to tumor formation. The theoretical approach is explicitly tested in the system with only two mutations using analytical calculations and Monte Carlo computer simulations. Excellent agreement with theoretical predictions is found for a large range of parameters, supporting our hypothesis and allowing us to better understand the mechanisms of cancer initiation. Our theoretical approach clarifies some important aspects of microscopic processes that lead to tumor formation.

Assessing, Testing and Estimating the Amount of Fine-Tuning by Means of Active Information

A general framework is introduced to estimate how much external information has been infused into a search algorithm, the so-called active information. This is rephrased as a test of fine-tuning, where tuning corresponds to the amount of pre-specified knowledge that the algorithm makes use of in order to reach a certain target. A function f quantifies specificity for each possible outcome x of a search, so that the target of the algorithm is a set of highly specified states, whereas fine-tuning occurs if it is much more likely for the algorithm to reach the target as intended than by chance. The distribution of a random outcome X of the algorithm involves a parameter θ that quantifies how much background information has been infused. A simple choice of this parameter is to use θf in order to exponentially tilt the distribution of the outcome of the search algorithm under the null distribution of no tuning, so that an exponential family of distributions is obtained. Such algorithms are obtained by iterating a Metropolis–Hastings type of Markov chain, which makes it possible to compute their active information under the equilibrium and non-equilibrium of the Markov chain, with or without stopping when the targeted set of fine-tuned states has been reached. Other choices of tuning parameters θ are discussed as well. Nonparametric and parametric estimators of active information and tests of fine-tuning are developed when repeated and independent outcomes of the algorithm are available. The theory is illustrated with examples from cosmology, student learning, reinforcement learning, a Moran type model of population genetics, and evolutionary programming.

The role of spatial structures of tissues in cancer initiation dynamics

It is widely believed that biological tissues evolved to lower the risks of cancer development. One of the specific ways to minimize the chances of tumor formation comes from proper spatial organization of tissues. However, the microscopic mechanisms of underlying processes remain not fully understood. We present a theoretical investigation on the role of spatial structures in cancer initiation dynamics. In our approach, the dynamics of single mutation fixations are analyzed using analytical calculations and computer simulations by mapping them to Moran processes on graphs with different connectivity that mimic various spatial structures. It is found that while the fixation probability is not affected by modifying the spatial structures of the tissues, the fixation times can change dramatically. The slowest dynamics is observed in "quasi-one-dimensional" structures, while the fastest dynamics is observed in "quasi-three-dimensional" structures. Theoretical calculations also suggest that there is a critical value of the degree of graph connectivity, which mimics the spatial dimension of the tissue structure, above which the spatial structure of the tissue has no effect on the mutation fixation dynamics. An effective discrete-state stochastic model of cancer initiation is utilized to explain our theoretical results and predictions. Our theoretical analysis clarifies some important aspects on the role of the tissue spatial structures in the cancer initiation processes.

Mathematical modeling for mutator phenotype and clonal selection advantage in the risk analysis of lung cancer

Cancer is one of the leading diseases for human mortality. Although substantial research works have been conducted to investigate the initiation and progression of cancer disease, it is still an active debate regarding the function of mutations conferring a clone advantage and the importance of mutator phenotypes caused by the mutation of stability genes. To address this issue further, we develop a mathematical model based on the incidence data of non-small cell lung cancer and small cell lung cancer from the Surveillance Epidemiology and End Results registry in the USA. The key biological parameters have been analyzed to investigate the potential effective measures for inhibiting the risk of lung cancer. Although the first event is the gene mutation that leads to clonal expansion of cells for lung cancer, the simulation results show that the clonal advantage of cancer cells alone is insufficient to cause tumorigenesis. Our analysis suggests that mutations in genes that keep genetic stability are critical in the development of lung cancer. This implies that mutator phenotype is an important indicator for the diagnosis of lung cancer, which can enable early detection and treatment to reduce the risk of lung cancer effectively. Furthermore, the parameter analysis indicates that it would be highly effective to control the risk of lung cancer by inhibiting the transformation rate from the normal cells to mutated cells and the clonal expansion of cells with fewer gene mutations.

Looking at Thyroid Cancer from the Tumor-Suppressor Genes Point of View

Thyroid cancer is the most frequent endocrine malignancy and accounts for approximately 1% of all diagnosed cancers. A variety of mechanisms are involved in the transformation of a normal tissue into a malignant one. Loss of tumor-suppressor gene (TSG) function is one of these mechanisms. The normal functions of TSGs include cell proliferation and differentiation control, genomic integrity maintenance, DNA damage repair, and signaling pathway regulation. TSGs are generally classified into three subclasses: (i) gatekeepers that encode proteins involved in cell cycle and apoptosis control; (ii) caretakers that produce proteins implicated in the genomic stability maintenance; and (iii) landscapers that, when mutated, create a suitable environment for malignant cell growth. Several possible mechanisms have been implicated in TSG inactivation. Reviewing the various TSG alteration types detected in thyroid cancers may help researchers to better understand the TSG defects implicated in the development/progression of this cancer type and to find potential targets for prognostic, predictive, diagnostic, and therapeutic purposes. Hence, the main purposes of this review article are to describe the various TSG inactivation mechanisms and alterations in human thyroid cancer, and the current therapeutic options for targeting TSGs in thyroid cancer.

Aspirin's effect on kinetic parameters of cells contributes to its role in reducing incidence of advanced colorectal adenomas, shown by a multiscale computational study

Aspirin intake has been shown to lead to significant protection against colorectal cancer, for example with an up to twofold reduction in colorectal adenoma incidence rates at higher doses. The mechanisms contributing to protection are not yet fully understood. While aspirin is an anti-inflammatory drug and can thus influence the tumor microenvironment, in vitro and in vivo experiments have recently shown that aspirin can also have a direct effect on cellular kinetics and fitness. It reduces the rate of tumor cell division and increases the rate of cell death. The question arises whether such changes in cellular fitness are sufficient to significantly contribute to the epidemiologically observed protection. To investigate this, we constructed a class of mathematical models of in vivo evolution of advanced adenomas, parameterized it with available estimates, and calculated population level incidence. Fitting the predictions to age incidence data revealed that only a model that included colonic crypt competition can account for the observed age-incidence curve. This model was then used to predict modified incidence patterns if cellular kinetics were altered as a result of aspirin treatment. We found that changes in cellular fitness that were within the experimentally observed ranges could reduce advanced adenoma incidence by a sufficient amount to account for age incidence data in aspirin-treated patient cohorts. While the mechanisms that contribute to the protective effect of aspirin are likely complex and multi-factorial, our study demonstrates that direct aspirin-induced changes of tumor cell fitness can significantly contribute to epidemiologically observed reduced incidence patterns.

OPTIMAL CONTROL OF GENETIC DIVERSITY IN THE MORAN MODEL WITH POPULATION GROWTH

In the Moran model of drift and selection of a mutant allele with population growth, instead of examining the consequences of pre-specified selection and population growth, the coexistence of the wild allele and the mutant allele becomes the maximization of the expected sojourn time in a given set. The process is controlled by the additional mortality of the mutant and by population growth. This makes it possible to retroactively assign fitness values as functions of the constraints, thus guiding a conservation policy or the achievement of a wishful proportion of mutants. This also gives the optimal conditions that have allowed an observed coexistence.

Emergent spatiotemporal population dynamics with cell-length control of synthetic microbial consortia

The increased complexity of synthetic microbial biocircuits highlights the need for distributed cell functionality due to concomitant increases in metabolic and regulatory burdens imposed on single-strain topologies. Distributed systems, however, introduce additional challenges since consortium composition and spatiotemporal dynamics of constituent strains must be robustly controlled to achieve desired circuit behaviors. Here, we address these challenges with a modeling-based investigation of emergent spatiotemporal population dynamics using cell-length control in monolayer, two-strain bacterial consortia. We demonstrate that with dynamic control of a strain's division length, nematic cell alignment in close-packed monolayers can be destabilized. We find that this destabilization confers an emergent, competitive advantage to smaller-length strains-but by mechanisms that differ depending on the spatial patterns of the population. We used complementary modeling approaches to elucidate underlying mechanisms: an agent-based model to simulate detailed mechanical and signaling interactions between the competing strains, and a reductive, stochastic lattice model to represent cell-cell interactions with a single rotational parameter. Our modeling suggests that spatial strain-fraction oscillations can be generated when cell-length control is coupled to quorum-sensing signaling in negative feedback topologies. Our research employs novel methods of population control and points the way to programming strain fraction dynamics in consortial synthetic biology.

Risk of lung cancer due to external environmental factor and epidemiological data analysis

Lung cancer is a cancer with the fastest growth in the incidence and mortality all over the world, which is an extremely serious threat to human's life and health. Evidences reveal that external environmental factors are the key drivers of lung cancer, such as smoking, radiation exposure and so on. Therefore, it is urgent to explain the mechanism of lung cancer risk due to external environmental factors experimentally and theoretically. However, it is still an open issue regarding how external environment factors affect lung cancer risk. In this paper, we summarize the main mathematical models involved the gene mutations for cancers, and review the application of the models to analyze the mechanism of lung cancer and the risk of lung cancer due to external environmental exposure. In addition, we apply the model described and the epidemiological data to analyze the influence of external environmental factors on lung cancer risk. The result indicates that radiation can cause significantly an increase in the mutation rate of cells, in particular the mutation in stability gene that leads to genomic instability. These studies not only can offer insights into the relationship between external environmental factors and human lung cancer risk, but also can provide theoretical guidance for the prevention and control of lung cancer.

Temporal order of mutations influences cancer initiation dynamics

Cancer is a set of genetic diseases that are driven by mutations. It was recently discovered that the temporal order of genetic mutations affects the cancer evolution and even the nature of the decease itself. The mechanistic origin of these observations, however, remain not well understood. Here we present a theoretical model for cancer initiation dynamics that allows us to quantify the impact of the temporal order of mutations. In our approach, the cancer initiation process is viewed as a set of stochastic transitions between discrete states defined by the different numbers of mutated cells. Using a first-passage analysis, probabilities and times before the cancer initiation are explicitly evaluated for two alternative sequences of two mutations. It is found that the probability of cancer initiation is determined only by the first mutation, while the dynamics depends on both mutations. In addition, it is shown that the acquisition of a mutation with higher fitness before mutation with lower fitness increases the probability of the tumor formation but delays the cancer initiation. Theoretical results are explained using effective free-energy landscapes.

Spectral analysis of transient amplifiers for death-birth updating constructed from regular graphs

A central question of evolutionary dynamics on graphs is whether or not a mutation introduced in a population of residents survives and eventually even spreads to the whole population, or becomes extinct. The outcome naturally depends on the fitness of the mutant and the rules by which mutants and residents may propagate on the network, but arguably the most determining factor is the network structure. Some structured networks are transient amplifiers. They increase for a certain fitness range the fixation probability of beneficial mutations as compared to a well-mixed population. We study a perturbation method for identifying transient amplifiers for death–birth updating. The method involves calculating the coalescence times of random walks on graphs and finding the vertex with the largest remeeting time. If the graph is perturbed by removing an edge from this vertex, there is a certain likelihood that the resulting perturbed graph is a transient amplifier. We test all pairwise nonisomorphic regular graphs up to a certain order and thus cover the whole structural range expressible by these graphs. For cubic and quartic regular graphs we find a sufficiently large number of transient amplifiers. For these networks we carry out a spectral analysis and show that the graphs from which transient amplifiers can be constructed share certain structural properties. Identifying spectral and structural properties may promote finding and designing such networks.

Mathematical modeling of multiple pathways in colorectal carcinogenesis using dynamical systems with Kronecker structure

Like many other types of cancer, colorectal cancer (CRC) develops through multiple pathways of carcinogenesis. This is also true for colorectal carcinogenesis in Lynch syndrome (LS), the most common inherited CRC syndrome. However, a comprehensive understanding of the distribution of these pathways of carcinogenesis, which allows for tailored clinical treatment and even prevention, is still lacking. We suggest a linear dynamical system modeling the evolution of different pathways of colorectal carcinogenesis based on the involved driver mutations. The model consists of different components accounting for independent and dependent mutational processes. We define the driver gene mutation graphs and combine them using the Cartesian graph product. This leads to matrix components built by the Kronecker sum and product of the adjacency matrices of the gene mutation graphs enabling a thorough mathematical analysis and medical interpretation. Using the Kronecker structure, we developed a mathematical model which we applied exemplarily to the three pathways of colorectal carcinogenesis in LS. Beside a pathogenic germline variant in one of the DNA mismatch repair (MMR) genes, driver mutations in APC, CTNNB1, KRAS and TP53 are considered. We exemplarily incorporate mutational dependencies, such as increased point mutation rates after MMR deficiency, and based on recent experimental data, biallelic somatic CTNNB1 mutations as common drivers of LS-associated CRCs. With the model and parameter choice, we obtained simulation results that are in concordance with clinical observations. These include the evolution of MMR-deficient crypts as early precursors in LS carcinogenesis and the influence of variants in MMR genes thereon. The proportions of MMR-deficient and MMR-proficient APC-inactivated crypts as first measure for the distribution among the pathways in LS-associated colorectal carcinogenesis are compatible with clinical observations. The approach provides a modular framework for modeling multiple pathways of carcinogenesis yielding promising results in concordance with clinical observations in LS CRCs.

On the waiting time until coordinated mutations get fixed in regulatory sequences

In this paper we consider the time evolution of a population of size N with overlapping generations, in the vicinity of m genes. We assume that this population is subject to point mutations, genetic drift, and selection. More specifically, we analyze the statistical distribution of the waiting time T_m until the expression of these genes have changed for all individuals, when transcription factors recognize and attach to short DNA-sequences (binding sites) within regulatory sequences in the neighborhoods of the m genes. The evolutionary dynamics is described by a multitype Moran process, where each individual is assigned a m×L regulatory array that consists of regulatory sequences with L nucleotides for all m genes. We study how the waiting time distribution depends on the number of genes, the mutation rate, the length of the binding sites, the length of the regulatory sequences, and the way in which the targeted binding sites are coordinated for different genes in terms of selection coefficients. These selection coefficients depend on how many binding sites have appeared so far, and possibly on their order of appearance. We also allow for back mutations, whereby some acquired binding sites may be lost over time. It is further assumed that the mutation rate is small enough to warrant a fixed state population, so that all individuals have the same regulatory array, at any given time point, until the next successful mutation arrives in some individual and spreads to the rest of the population. We further incorporate stochastic tunneling, whereby successful mutations get mutated before their fixation. A crucial part of our approach is to divide the huge state space of regulatory arrays into a small number of components, assuming that the array component varies as a Markov process over time. This implies that T_m is the time until this Markov process hits an absorbing state, with a phase-type distribution. A number of interesting results can be derived from our general setup, for instance that the expected waiting time increases exponentially with m, for a selectively neutral model, when back-mutations are possible.

Engineering in Medicine To Address the Challenge of Cancer Drug Resistance: From Micro- and Nanotechnologies to Computational and Mathematical Modeling

Drug resistance has profoundly limited the success of cancer treatment, driving relapse, metastasis, and mortality. Nearly all anticancer drugs and even novel immunotherapies, which recalibrate the immune system for tumor recognition and destruction, have succumbed to resistance development. Engineers have emerged across mechanical, physical, chemical, mathematical, and biological disciplines to address the challenge of drug resistance using a combination of interdisciplinary tools and skill sets. This review explores the developing, complex, and under-recognized role of engineering in medicine to address the multitude of challenges in cancer drug resistance. Looking through the "lens" of intrinsic, extrinsic, and drug-induced resistance (also referred to as "tolerance"), we will discuss three specific areas where active innovation is driving novel treatment paradigms: (1) nanotechnology, which has revolutionized drug delivery in desmoplastic tissues, harnessing physiochemical characteristics to destroy tumors through photothermal therapy and rationally designed nanostructures to circumvent cancer immunotherapy failures, (2) bioengineered tumor models, which have benefitted from microfluidics and mechanical engineering, creating a paradigm shift in physiologically relevant environments to predict clinical refractoriness and enabling platforms for screening drug combinations to thwart resistance at the individual patient level, and (3) computational and mathematical modeling, which blends in silico simulations with molecular and evolutionary principles to map mutational patterns and model interactions between cells that promote resistance. On the basis that engineering in medicine has resulted in discoveries in resistance biology and successfully translated to clinical strategies that improve outcomes, we suggest the proliferation of multidisciplinary science that embraces engineering.

Using statistical methods to model the fine-tuning of molecular machines and systems

Fine-tuning has received much attention in physics, and it states that the fundamental constants of physics are finely tuned to precise values for a rich chemistry and life permittance. It has not yet been applied in a broad manner to molecular biology. However, in this paper we argue that biological systems present fine-tuning at different levels, e.g. functional proteins, complex biochemical machines in living cells, and cellular networks. This paper describes molecular fine-tuning, how it can be used in biology, and how it challenges conventional Darwinian thinking. We also discuss the statistical methods underpinning fine-tuning and present a framework for such analysis.

A Survey of Evolutionary Games in Biology

The applications of game theory in biology are numerous and include the comparison and modeling situations between bacteria, viruses, etc. This work provides insights about the connection between biology and evolving populations with classical and quantum evolutionary game theory and explains the benefits of unconventional computing methods in the study of such phenomena. In particular, the introduction of automata brings new possibilities into the decision-making process.

Mathematical modelling the pathway of genomic instability in lung cancer

Genomic instability plays a significant role in lung cancer. Although substantial research has been conducted using both clinical and theoretical studies, it is still a hotly debated issue to whether genomic instability is necessary or whether genomic instability precedes oncogenes activation and tumor suppressor genes inactivation for lung cancer. In response to this issue, we come up with a mathematical model incorporating effects of genomic instability to investigate the genomic instability pathway of human lung cancer. The presented model are applied to match the incidence rate data of lung cancer from the Life Span Study cohort of the atomic bomb survivors in Nagasaki and Hiroshima and the Surveillance Epidemiology and End Results registry in the United States. Model results suggest that genomic instability is necessary in the tumorigenesis of lung cancer, and genomic instability has no significant impact on the net proliferation rate of cells by statistical criteria. By comparing the results of the LSS data to those of the SEER data, we conclude that the genomic instability pathway exhibits a sensitivity to radiation exposure, more intensive in male patients.

Passenger mutations can accelerate tumour suppressor gene inactivation in cancer evolution

Carcinogenesis is an evolutionary process whereby cells accumulate multiple mutations. Besides the 'driver mutations' that cause the disease, cells also accumulate a number of other mutations with seemingly no direct role in this evolutionary process. They are called passenger mutations. While it has been argued that passenger mutations render tumours more fragile due to reduced fitness, the role of passenger mutations remains understudied. Using evolutionary computational models, we demonstrate that in the context of tumour suppressor gene inactivation (and hence fitness valley crossing), the presence of passenger mutations can accelerate the rate of evolution by reducing overall population fitness and increasing the relative fitness of intermediate mutants in the fitness valley crossing pathway. Hence, the baseline rate of tumour suppressor gene inactivation might be faster than previously thought. Conceptually, parallels are found in the field of turbulence and pattern formation, where instabilities can be driven by perturbations that are damped (disadvantageous), but provide a richer set of pathways such that a system can achieve some desired goal more readily. This highlights, through a number of novel parallels, the relevance of physical sciences in oncology.

The expected time to cross extended fitness plateaus

For a population to acquire a complex adaptation requiring multiple individually neutral mutations, it must cross a plateau in the fitness landscape. We consider plateaus involving three mutations, and show that large populations can cross them rapidly via lineages that acquire multiple mutations while remaining at low frequency, much faster than the ∝μ³ rate for simultaneous triple mutations. Plateau-crossing is fastest for very large populations. At intermediate population sizes, recombination can greatly accelerate adaptation by combining independent mutant lineages to form triple-mutants. For more frequent recombination, such that the population is kept near linkage equilibrium, we extend our analysis to find simple expressions for the expected time to cross plateaus of arbitrary width.

Multiple infection of cells changes the dynamics of basic viral evolutionary processes

The infection of cells by multiple copies of a given virus can impact viral evolution in a variety of ways, yet some of the most basic evolutionary dynamics remain underexplored. Using computational models, we investigate how infection multiplicity affects the fixation probability of mutants, the rate of mutant generation, and the timing of mutant invasion. An important insight from these models is that for neutral and disadvantageous phenotypes, rare mutants initially enjoy a fitness advantage in the presence of multiple infection of cells. This arises because multiple infection allows the rare mutant to enter more target cells and to spread faster, while it does not accelerate the spread of the resident wild-type virus. The rare mutant population can increase by entry into both uninfected and wild-type-infected cells, while the established wild-type population can initially only grow through entry into uninfected cells. Following this initial advantageous phase, the dynamics are governed by drift or negative selection, respectively, and a higher multiplicity reduces the chances that mutants fix in the population. Hence, while increased infection multiplicity promotes the presence of neutral and disadvantageous mutants in the short-term, it makes it less likely in the longer term. We show how these theoretical insights can be useful for the interpretation of experimental data on virus evolution at low and high multiplicities. The dynamics explored here provide a basis for the investigation of more complex viral evolutionary processes, including recombination, reassortment, as well as complementary/inhibitory interactions.

The role of telomere shortening in carcinogenesis: A hybrid stochastic-deterministic approach

Genome instability is a characteristic of most cancers, contributing to the acquisition of genetic alterations that drive tumor progression. One important source of genome instability is linked to telomere dysfunction in cells with critically short telomeres that lack p53-mediated surveillance of genomic integrity. Here we research the probability that cancer emerges through an evolutionary pathway that includes a telomere-induced phase of genome instability. To implement our models we use a hybrid stochastic-deterministic approach, which allows us to perform large numbers of simulations using biologically realistic population sizes and mutation rates, circumventing the traditional limitations of fully stochastic algorithms. The hybrid methodology should be easily adaptable to a wide range of evolutionary problems. In particular, we model telomere shortening and the acquisition of two mutations: Telomerase activation and p53 inactivation. We find that the death rate of unstable cells, and the number of cell divisions that p53 mutants can sustain beyond the normal senescence setpoint determine the likelihood that the first double mutant originates in a cell with telomere-induced instability. The model has applications to an influential telomerase-null mouse model and p16 silenced human cells. We end by discussing algorithmic performance and a measure for the accuracy of the hybrid approximation.

Fitness variation in isogenic populations leads to a novel evolutionary mechanism for crossing fitness valleys

Individuals in a population often have different fitnesses even when they have identical genotypes, but the effect of this variation on the evolution of a population through complicated fitness landscapes is unknown. Here, we investigate how populations with non-genetic fitness variation cross fitness valleys, common barriers to adaptation in rugged fitness landscapes in which a population must pass through a deleterious intermediate to arrive at a final advantageous stage. We develop a stochastic computational model describing the dynamics of an asexually reproducing population crossing a fitness valley, in which individuals of the same evolutionary stage can have variable fitnesses. We find that fitness variation that persists over multiple generations increases the rate of valley crossing through a novel evolutionary mechanism different from previously characterized mechanisms such as stochastic tunneling. By reducing the strength of selection against deleterious intermediates, persistent fitness variation allows for faster adaptation through rugged fitness landscapes.

Key issues review: evolution on rugged adaptive landscapes

Adaptive landscapes represent a mapping between genotype and fitness. Rugged adaptive landscapes contain two or more adaptive peaks: allele combinations with higher fitness than any of their neighbors in the genetic space. How do populations evolve on such rugged landscapes? Evolutionary biologists have struggled with this question since it was first introduced in the 1930s by Sewall Wright. Discoveries in the fields of genetics and biochemistry inspired various mathematical models of adaptive landscapes. The development of landscape models led to numerous theoretical studies analyzing evolution on rugged landscapes under different biological conditions. The large body of theoretical work suggests that adaptive landscapes are major determinants of the progress and outcome of evolutionary processes. Recent technological advances in molecular biology and microbiology allow experimenters to measure adaptive values of large sets of allele combinations and construct empirical adaptive landscapes for the first time. Such empirical landscapes have already been generated in bacteria, yeast, viruses, and fungi, and are contributing to new insights about evolution on adaptive landscapes. In this Key Issues Review we will: (i) introduce the concept of adaptive landscapes; (ii) review the major theoretical studies of evolution on rugged landscapes; (iii) review some of the recently obtained empirical adaptive landscapes; (iv) discuss recent mathematical and statistical analyses motivated by empirical adaptive landscapes, as well as provide the reader with instructions and source code to implement simulations of evolution on adaptive landscapes; and (v) discuss possible future directions for this exciting field.

The role of backward cell migration in two-hit mutants' production in the stem cell niche

It has been discovered that there are two stem cell groups in the intestinal crypts: central stem cells (CeSCs), which are at the very bottom of the crypt, and border stem cells (BSCs), which are located between CeSCs and transit amplifying cells (TAs). Moreover, backward cell migration from BSCs to CeSCs has been observed. Recently, a bi-compartmental stochastic model, which includes CeSCs and BSCs, has been developed to investigate the probability of two-hit mutant production in the stem cell niche. In this project, we improve this stochastic model by adding the probability of backward cell migration to the model. The model suggests that the probability of two-hit mutant production increases when the frequency of backward cell migration increases. Furthermore, a small non-zero probability of backward cell migration leads to the largest range of optimal values for the frequency of symmetric divisions and the portion of divisions at each stem cell compartment in terms of delaying 2-hit mutant production. Moreover, the probability of two-hit mutant production is more sensitive to the probability of symmetric divisions than to the rate of backward cell migrations. The highest probability of two-hit mutant production corresponds to the case when all stem cell’s divisions are asymmetric.

Effect of Population Size and Mutation Rate on the Evolution of RNA Sequences on an Adaptive Landscape Determined by RNA Folding

The dynamics of populations evolving on an adaptive landscape depends on multiple factors, including the structure of the landscape, the rate of mutations, and effective population size. Existing theoretical work often makes ad hoc and simplifying assumptions about landscape structure, whereas experimental work can vary important parameters only to a limited extent. We here overcome some of these limitations by simulating the adaptive evolution of RNA molecules, whose fitness is determined by the thermodynamics of RNA secondary structure folding. We study the influence of mutation rates and population sizes on final mean population fitness, on the substitution rates of mutations, and on population diversity. We show that evolutionary dynamics cannot be understood as a function of mutation rate µ, population size N, or population mutation rate Nµ alone. For example, at a given mutation rate, clonal interference prevents the fixation of beneficial mutations as population size increases, but larger populations still arrive at a higher mean fitness. In addition, at the highest population mutation rates we study, mean final fitness increases with population size, because small populations are driven to low fitness by the relatively higher incidence of mutations they experience. Our observations show that mutation rate and population size can interact in complex ways to influence the adaptive dynamics of a population on a biophysically motivated fitness landscape.

Forward and backward evolutionary processes and allele frequency spectrum in a cancer cell population

A cancer grows from a single cell, thereby constituting a large cell population. In this work, we are interested in how mutations accumulate in a cancer cell population. We provide a theoretical framework of the stochastic process in a cancer cell population and obtain near exact expressions of allele frequency spectrum or AFS (only continuous approximation is involved) from both forward and backward treatments under a simple setting; all cells undergo cell divisions and die at constant rates, b and d, respectively, such that the entire population grows exponentially. This setting means that once a parental cancer cell is established, in the following growth phase, all mutations are assumed to have no effect on b or d (i.e., neutral or passengers). Our theoretical results show that the difference from organismal population genetics is mainly in the coalescent time scale, and the mutation rate is defined per cell division, not per time unit (e.g., generation). Except for these two factors, the basic logic is very similar between organismal and cancer population genetics, indicating that a number of well established theories of organismal population genetics could be translated to cancer population genetics with simple modifications.

Suppressors of selection

Inspired by recent works on evolutionary graph theory, an area of growing interest in mathematical and computational biology, we present examples of undirected structures acting as suppressors of selection for any fitness value r > 1. This means that the average fixation probability of an advantageous mutant or invader individual placed at some node is strictly less than that of this individual placed in a well-mixed population. This leads the way to study more robust structures less prone to invasion, contrary to what happens with the amplifiers of selection where the fixation probability is increased on average for advantageous invader individuals. A few families of amplifiers are known, although some effort was required to prove it. Here, we use computer aided techniques to find an exact analytical expression of the fixation probability for some graphs of small order (equal to 6, 8 and 10) proving that selection is effectively reduced for r > 1. Some numerical experiments using Monte Carlo methods are also performed for larger graphs and some variants.

With a little help from my friends: cooperation can accelerate the rate of adaptive valley crossing

Background

Natural selection favors changes that lead to genotypes possessing high fitness. A conflict arises when several mutations are required for adaptation, but each mutation is separately deleterious. The process of a population evolving from a genotype encoding for a local fitness maximum to a higher fitness genotype is termed an adaptive peak shift.

Results

Here we suggest cooperative behavior as a factor that can facilitate adaptive peak shifts. We model cooperation in a public goods scenario, wherein each individual contributes resources that are later equally redistributed among all cooperating individuals. We use mathematical modeling and stochastic simulations to study the effect of cooperation on peak shifts in both panmictic and structured populations. Our results show that cooperation can substantially affect the rate of complex adaptation. Furthermore, we show that cooperation increases the population diversity throughout the peak shift process, thus increasing the robustness of the population to sudden environmental changes.

Conclusions

We provide a new explanation to adaptive valley crossing in natural populations and suggest that the long term evolution of a species depends on its social behavior.

Electronic supplementary material

The online version of this article (doi:10.1186/s12862-017-0983-2) contains supplementary material, which is available to authorized users.

CellTrans: An R Package to Quantify Stochastic Cell State Transitions

Many normal and cancerous cell lines exhibit a stable composition of cells in distinct states which can, e.g., be defined on the basis of cell surface markers. There is evidence that such an equilibrium is associated with stochastic transitions between distinct states. Quantifying these transitions has the potential to better understand cell lineage compositions. We introduce CellTrans, an R package to quantify stochastic cell state transitions from cell state proportion data from fluorescence-activated cell sorting and flow cytometry experiments. The R package is based on a mathematical model in which cell state alterations occur due to stochastic transitions between distinct cell states whose rates only depend on the current state of a cell. CellTrans is an automated tool for estimating the underlying transition probabilities from appropriately prepared data. We point out potential analytical challenges in the quantification of these cell transitions and explain how CellTrans handles them. The applicability of CellTrans is demonstrated on publicly available data on the evolution of cell state compositions in cancer cell lines. We show that CellTrans can be used to (1) infer the transition probabilities between different cell states, (2) predict cell line compositions at a certain time, (3) predict equilibrium cell state compositions, and (4) estimate the time needed to reach this equilibrium. We provide an implementation of CellTrans in R, freely available via GitHub (https://github.com/tbuder/CellTrans).

Genetic load makes cancer cells more sensitive to common drugs: evidence from Cancer Cell Line Encyclopedia

Genetic alterations initiate tumors and enable the evolution of drug resistance. The pro-cancer view of mutations is however incomplete, and several studies show that mutational load can reduce tumor fitness. Given its negative effect, genetic load should make tumors more sensitive to anticancer drugs. Here, we test this hypothesis across all major types of cancer from the Cancer Cell Line Encyclopedia, which provides genetic and expression data of 496 cell lines together with their response to 24 common anticancer drugs. We found that the efficacy of 9 out of 24 drugs showed significant association with genetic load in a pan-cancer analysis. The associations for some tissue-drug combinations were remarkably strong, with genetic load explaining up to 83% of the variance in the drug response. Overall, the role of genetic load depended on both the drug and the tissue type with 10 tissues being particularly vulnerable to genetic load. We also identified changes in gene expression associated with increased genetic load, which included cell-cycle checkpoints, DNA damage and apoptosis. Our results show that genetic load is an important component of tumor fitness and can predict drug sensitivity. Beyond being a biomarker, genetic load might be a new, unexplored vulnerability of cancer.

Modeling dynamics of mutants in heterogeneous stem cell niche

Studying the stem cell (SC) niche architecture is a crucial step for investigating the process of oncogenesis and obtaining an effective stem cell therapy for various cancers. Recently, it has been observed that there are two groups of SCs in the SC niche collaborating with each other to maintain tissue homeostasis: border stem cells (BSCs), which are responsible in controlling the number of non-stem cells as well as stem cells, and central stem cells (CeSCs), which regulate the SC niche. Here, we develop a bi-compartmental stochastic model for the SC niche to study the spread of mutants within the niche. The analytic calculations and numeric simulations, which are in perfect agreement, reveal that in order to delay the spread of mutants in the SC niche, a small but non-zero number of SC proliferations must occur in the CeSC compartment. Moreover, the migration of BSCs to CeSCs delays the spread of mutants. Furthermore, the fixation probability of mutants in the SC niche is independent of types of SC division as long as all SCs do not divide fully asymmetrically. Additionally, the progeny of CeSCs have a much higher chance than the progeny of BSCs to take over the entire niche.

A Hierarchical Probability Model of Colon Cancer

We consider a model of fixed sizeN= 2lin which there arelgenerations of daughter cells and a stem cell. In each generationithere are 2i−1daughter cells. At each integral time unit the cells split so that the stem cell splits into a stem cell and generation 1 daughter cell and the generationidaughter cells become two cells of generationi+1. The last generation is removed from the population. A stem cell acquires first and second mutations at ratesu1andu2, and a daughter cell acquires first and second mutations at ratesv1andv2. We find the distribution for the time it takes to acquire two mutations asNgoes to ∞ and the mutation rates go to 0. The mutation rates may tend to 0 at different speeds. We also find the distribution for the locations of the mutations. In particular, we determine whether or not the mutations occur on a stem cell and if not, at what generation in the daughter cells they occur. Several outcomes are possible, depending on how fast the rates go to 0. The model considered has been proposed by Komarova (2007) as a model for colon cancer.

Replicator dynamics of cancer stem cell: Selection in the presence of differentiation and plasticity

The cancer stem cell hypothesis has evolved into one of the most important paradigms in cancer research. According to cancer stem cell hypothesis, somatic mutations in a subpopulation of cells can transform them into cancer stem cells with the unique potential of tumour initiation. Stem cells have the potential to produce lineages of non-stem cell populations (differentiated cells) via a ubiquitous hierarchal division scheme. Differentiation of a stem cell into (partially) differentiated cells can happen either symmetrically or asymmetrically. The selection dynamics of a mutant cancer stem cell should be investigated in the light of a stem cell proliferation hierarchy and presence of a non-stem cell population. By constructing a three-compartment Moran-type model composed of normal stem cells, mutant (cancer) stem cells and differentiated cells, we derive the replicator dynamics of stem cell frequencies where asymmetric differentiation and differentiated cell death rates are included in the model. We determine how these new factors change the conditions for a successful mutant invasion and discuss the variation on the steady state fraction of the population as different model parameters are changed. By including the phenotypic plasticity/dedifferentiation, in which a progenitor/differentiated cell can transform back into a cancer stem cell, we show that the effective fitness of mutant stem cells is not only determined by their proliferation and death rates but also according to their dedifferentiation potential. By numerically solving the model we derive the phase diagram of the advantageous and disadvantageous phases of cancer stem cells in the space of proliferation and dedifferentiation potentials. The result shows that at high enough dedifferentiation rates even a previously disadvantageous mutant can take over the population of normal stem cells. This observation has implications in different areas of cancer research including experimental observations that imply metastatic cancer stem cell types might have lower proliferation potential than other stem cell phenotypes while showing much more phenotypic plasticity and can undergo clonal expansion.

Model-Based Evaluation of Spontaneous Tumor Regression in Pilocytic Astrocytoma

Pilocytic astrocytoma (PA) is the most common brain tumor in children. This tumor is usually benign and has a good prognosis. Total resection is the treatment of choice and will cure the majority of patients. However, often only partial resection is possible due to the location of the tumor. In that case, spontaneous regression, regrowth, or progression to a more aggressive form have been observed. The dependency between the residual tumor size and spontaneous regression is not understood yet. Therefore, the prognosis is largely unpredictable and there is controversy regarding the management of patients for whom complete resection cannot be achieved. Strategies span from pure observation (wait and see) to combinations of surgery, adjuvant chemotherapy, and radiotherapy. Here, we introduce a mathematical model to investigate the growth and progression behavior of PA. In particular, we propose a Markov chain model incorporating cell proliferation and death as well as mutations. Our model analysis shows that the tumor behavior after partial resection is essentially determined by a risk coefficient γ, which can be deduced from epidemiological data about PA. Our results quantitatively predict the regression probability of a partially resected benign PA given the residual tumor size and lead to the hypothesis that this dependency is linear, implying that removing any amount of tumor mass will improve prognosis. This finding stands in contrast to diffuse malignant glioma where an extent of resection threshold has been experimentally shown, below which no benefit for survival is expected. These results have important implications for future therapeutic studies in PA that should include residual tumor volume as a prognostic factor.

Modeling Invasion Dynamics with Spatial Random-Fitness Due to Micro-Environment

Numerous experimental studies have demonstrated that the microenvironment is a key regulator influencing the proliferative and migrative potentials of species. Spatial and temporal disturbances lead to adverse and hazardous microenvironments for cellular systems that is reflected in the phenotypic heterogeneity within the system. In this paper, we study the effect of microenvironment on the invasive capability of species, or mutants, on structured grids (in particular, square lattices) under the influence of site-dependent random proliferation in addition to a migration potential. We discuss both continuous and discrete fitness distributions. Our results suggest that the invasion probability is negatively correlated with the variance of fitness distribution of mutants (for both advantageous and neutral mutants) in the absence of migration of both types of cells. A similar behaviour is observed even in the presence of a random fitness distribution of host cells in the system with neutral fitness rate. In the case of a bimodal distribution, we observe zero invasion probability until the system reaches a (specific) proportion of advantageous phenotypes. Also, we find that the migrative potential amplifies the invasion probability as the variance of fitness of mutants increases in the system, which is the exact opposite in the absence of migration. Our computational framework captures the harsh microenvironmental conditions through quenched random fitness distributions and migration of cells, and our analysis shows that they play an important role in the invasion dynamics of several biological systems such as bacterial micro-habitats, epithelial dysplasia, and metastasis. We believe that our results may lead to more experimental studies, which can in turn provide further insights into the role and impact of heterogeneous environments on invasion dynamics.

Can Peto's paradox be used as the null hypothesis to identify the role of evolution in natural resistance to cancer? A critical review

BACKGROUND

Carcinogenesis affects not only humans but almost all metazoan species. Understanding the rules driving the occurrence of cancers in the wild is currently expected to provide crucial insights into identifying how some species may have evolved efficient cancer resistance mechanisms. Recently the absence of correlation across species between cancer prevalence and body size (coined as Peto's paradox) has attracted a lot of attention. Indeed, the disparity between this null hypothesis, where every cell is assumed to have an identical probability to undergo malignant transformation, and empirical observations is particularly important to understand, due to the fact that it could facilitate the identification of animal species that are more resistant to carcinogenesis than expected. Moreover it would open up ways to identify the selective pressures that may be involved in cancer resistance. However, Peto's paradox relies on several questionable assumptions, complicating the interpretation of the divergence between expected and observed cancer incidences.

DISCUSSIONS

Here we review and challenge the different hypotheses on which this paradox relies on with the aim of identifying how this null hypothesis could be better estimated in order to provide a standard protocol to study the deviation between theoretical/theoretically predicted and observed cancer incidence. We show that due to the disproportion and restricted nature of available data on animal cancers, applying Peto's hypotheses at species level could result in erroneous conclusions, and actually assume the existence of a paradox. Instead of using species level comparisons, we propose an organ level approach to be a more accurate test of Peto's assumptions.

SUMMARY

The accuracy of Peto's paradox assumptions are rarely valid and/or quantifiable, suggesting the need to reconsider the use of Peto's paradox as a null hypothesis in identifying the influence of natural selection on cancer resistance mechanisms.

Steering Evolution with Sequential Therapy to Prevent the Emergence of Bacterial Antibiotic Resistance

The increasing rate of antibiotic resistance and slowing discovery of novel antibiotic treatments presents a growing threat to public health. Here, we consider a simple model of evolution in asexually reproducing populations which considers adaptation as a biased random walk on a fitness landscape. This model associates the global properties of the fitness landscape with the algebraic properties of a Markov chain transition matrix and allows us to derive general results on the non-commutativity and irreversibility of natural selection as well as antibiotic cycling strategies. Using this formalism, we analyze 15 empirical fitness landscapes of E. coli under selection by different β-lactam antibiotics and demonstrate that the emergence of resistance to a given antibiotic can be either hindered or promoted by different sequences of drug application. Specifically, we demonstrate that the majority, approximately 70%, of sequential drug treatments with 2-4 drugs promote resistance to the final antibiotic. Further, we derive optimal drug application sequences with which we can probabilistically 'steer' the population through genotype space to avoid the emergence of resistance. This suggests a new strategy in the war against antibiotic-resistant organisms: drug sequencing to shepherd evolution through genotype space to states from which resistance cannot emerge and by which to maximize the chance of successful therapy.