Modeling Evolution of Heterogeneous Populations: Theory and Applications

By Irina Kareva and Georgy Karev

Modeling Evolution of Heterogeneous Populations: Theory and Applications describes, develops and provides applications of a method that allows incorporating population heterogeneity into systems of ordinary and discrete differential equations without significantly increasing system dimensionality. The method additionally allows making use of results of bifurcation analysis performed on simplified homogeneous systems, thereby building on the existing body of tools and knowledge and expanding applicability and predictive power of many mathematical models.

• Introduces Hidden Keystone Variable (HKV) method, which allows modeling evolution of heterogenous populations, while reducing multi-dimensional selection systems to low-dimensional systems of differential equations
• Demonstrates that replicator dynamics is governed by the principle of maximal relative entropy that can be derived from the dynamics of selection systems instead of being postulated
• Discusses mechanisms behind models of both Darwinian and non-Darwinian selection
• Provides examples of applications to various fields, including cancer growth, global demography, population extinction, tragedy of the commons and resource sustainability, among others
• Helps inform differences in underlying mechanisms of population growth from experimental observations, taking one from experiment to theory and back

• Language: English
• Publisher: Academic Press
• Release date: Oct 16, 2019
• ISBN: 9780128144329

Irina Kareva

Dr. Irina Kareva is a theoretical biologist, and the primary focus of her research involves using mathematical modeling to study cancer as an evolving ecosystem within the human body, where heterogeneous populations of cancer cells compete for limited resources (i.e., oxygen and glucose), cooperate with each other to fight off predators (the immune system), and disperse and migrate (metastases). In 2017 Dr. Kareva gave a TED talk on using mathematical modeling for biological research. Dr. Kareva's book Understanding cancer from a systems biology point of view: from observation to theory and back was published by Elsevier in 2018. Dr. Kareva is a Senior Scientist in Simulation and Modeling at EMD Serono, Merck KGaA, where she develops quantitative systems pharmacology (QSP) models to help understand and predict dynamics of new therapeutics.

Book preview

mathematician

Chapter 1

Using mathematical modeling to ask meaningful biological questions

Abstract

Classical approaches to analyzing dynamical systems, including bifurcation analysis, can provide invaluable insights into the underlying structure of a mathematical model and the spectrum of all possible dynamical behaviors. However, these models frequently fail to take into account population heterogeneity. While heterogeneity is critically important to understanding and predicting the behavior of any evolving system, this characteristic is commonly omitted when analyzing many mathematical models of ecological systems. Attempts to include population heterogeneity frequently result in expanding system dimensionality, effectively preventing qualitative analysis. However, Reduction theorem, or hidden keystone variable (HKV) method, allows incorporating population heterogeneity while still permitting the use of classical bifurcation analysis. A combination of these methods allows visualization of evolutionary trajectories and permits making meaningful predictions about dynamics over time of evolving populations.

Keywords

HKV method; Population heterogeneity; Replicator equations; Selection system

1.1 Introduction

Heterogeneity is a major driving force behind the dynamics of evolving systems. When it is heritable and when it affects fitness, heterogeneity is what makes evolution possible (Bell, 2008; Darwin, 1880; Johnson, 1976; Page, 2010). This comes from the fact that the environment in which the individuals interact is composed not only of the outside world (such as the resources necessary for survival, or members of other species) but also of individuals themselves. Therefore, selective pressures that are imposed on the individuals come both from the environment and from one another. Furthermore, selective pressures that individuals experience from one another will be imposed and perceived differently depending on population composition, which in turn may be changing as a result of these selective pressures. This selection process underlies Dobzhansky's famous thesis: nothing in biology makes sense except in light of evolution (Dobzhansky, 1973).

In a majority of conceptual and often even descriptive mathematical models of population dynamics—whether it be models of predator-prey interactions, spread of infectious diseases or tumor growth—population homogeneity is the first simplification that is made. The population is not treated as homogeneous per se; rather, one assumes an average rate of growth, death, or infectiousness as a reasonable enough approximation if the system has reached a quasi-stable state of evolutionary development. However, by ignoring population heterogeneity in such a way, one ends up either ignoring natural selection or assuming that it has already done its work. This assumption is often incorrect within the context of such models however, since natural selection may in fact be a key driver behind the dynamics of systems that are of most interest and importance.

The process of evolution of heterogeneous populations is typically modeled using replicator equations that capture the basic tenet of Darwinism (Hofbauer and Sigmund, 1998; Nowak and Sigmund, 2004). These equations capture the process of selection in heterogeneous populations as they demonstrate how when an individual's value of a heritable characteristic is above average, that individual stays in the population, but when the heritable characteristic is below average, the individual is expunged. However, such models have a major drawback: modeling high levels of heterogeneity is accompanied by an inevitable increase of system dimensionality, which makes obtaining any kind of qualitative understanding of the system nearly impossible. Assuming population homogeneity makes systems of equations computationally and sometimes even analytically manageable, but at the cost of losing many of the system dynamics caused by intraspecies interactions and natural selection.

Despite their shortcomings, parametrically homogeneous systems can still provide exceptionally valuable information about the structure of the system through the use of extensively developed analytical techniques, such as bifurcation analysis (Kuznetsov, 2013). A skillfully constructed bifurcation diagram can reveal various possible dynamical regimes of a system that result from variations in parameter values and initial conditions, and also provide analytical boundaries as functions of system parameters. This information can then be used to construct a theoretical framework for understanding a biological system that could never have been obtained experimentally.

In this book, we will describe in detail a method to introduce population heterogeneity back into equation-based models using the Reduction theorem. Also known as parameter distribution technique or hidden keystone variable (HKV) method, Reduction theorem can use and build on insights obtained from bifurcation analysis, while incorporating population heterogeneity. Reduction theorem allows the dynamics of an evolving system to be investigated more fully, while overcoming the problem of immense system dimensionality in a wide class of mathematical models.

1.2 General strategy

The key steps for using the HKV method are as follows. Assume a population of individuals is composed at time t of clones x(t, a). Each individual clone x(t, a) is characterized by parameter value a ∈ Aif the system is discrete, and N(t) = ∫Ax(t, a)da can change over time due to system dynamics. Consequently the mean value of the parameter Et[a] now becomes a function of time and changes over time as well.

In the upcoming chapters we show how to analyze a parametrically heterogeneous system using the following steps:

1.Analyze the autonomous parametrically homogeneous system to the extent possible using well-developed analytical tools, such as bifurcation analysis.

2.Replace parameter a with its mean value Et[a], which is a function of time.

3.Introduce an auxiliary system of differential equations to define keystone variables that determine the actual dynamics of the system. (Note: the term keystone is used here to parallel the function of keystone species in ecology. Just like keystone species have disproportionately large effect on their environment relative to their abundance, keystone variables determine the direction in which the system will evolve, without being explicitly present in the original system.)

4.Express the distribution of the distributed parameter through keystone variables. This transformation allows finding all statistical characteristics of interest, including parameter's mean and variance, which now change over time due to system dynamics. The mean of the parameter can now travel through the different domains of the phase-parameter portrait of the original parametrically homogeneous system.

5.Calculate numerical solutions.

Exact formulation of the Reduction theorem and the theory underlying the HKV method will be found later in the book. A summary definitions and associated notation are provided in Table 1.

Table 1

1.3 Advantages and drawbacks of the Reduction theorem

One of the most important properties of this method is that it allows reducing an otherwise many- or even infinitely-dimensional system to low dimensionality.

However, as with any method, this method is not universal. Most importantly, the transformation can be done (with some generalizations) only to Kolmogorov type equations of the form x(t)'= x(t)F(t, Et[f(a)]), where:

•x(t) is a vector,

•a is a parameter or a vector of parameters that characterize individual heterogeneity within the population,

•Et[f(a)] is of system-specific form.

Reduction theorem can also increase the dimensionality of the original parametrically homogeneous system at a possible cost of auxiliary keystone equations (although these would typically be up to only one or two extra equations, depending on the original system). Finally, the resulting system is typically non-autonomous, so one cannot perform standard bifurcation analysis.

When studying numerical solutions of such parametrically heterogeneous systems, trajectories can be observed that could not previously have been seen in parametrically homogeneous systems. This phenomenon arises from the expected value of the parameter traveling through the phase parameter portrait, and the system undergoes corresponding qualitative phase transition as the parameter's expected value crosses the bifurcation boundaries. Furthermore, if there exists a complete bifurcation diagram for the specific parametrically homogeneous model, the boundaries crossed during system evolution can be identified analytically.

Classic techniques for analyzing dynamical systems, such as bifurcation theory (Kuznetsov, 2013), can provide critical insights into the possible dynamical regimes that a system can realize. Unfortunately, doing full bifurcation analysis is labor intensive and is not always possible due to increasing complexities of many models. However, a very rich body of literature exists of fully analyzed parametrically homogeneous models in many fields, including ecology (Bazykin, 1998; Berezovskaya et al., 2005), epidemiology (Brauer and Castillo-Chavez, 2001), among others. As the examples presented throughout this book will demonstrate, even relatively simple two-dimensional systems can reveal rich, unexpected and meaningful behaviors. Application of the HKV-method to introduce population heterogeneity in a meaningful way and utilizing previously performed analysis can reveal a new layer of understanding of many existing models that was not accessible before.

In Chapter 1 of this book, we introduce the HKV method for modeling population heterogeneity. Chapter 2 shows how applications of these methods to some classical models can reveal new and unexpected dynamical behaviors. Chapter 3 demonstrates how the HKV-method can be applied to more complex biological systems, including models of world demography, microbial resistance to antibiotics, and the dynamics of tree stand self-thinning. In Chapter 4, we go over more in-depth theory. Chapter 5 features numerous examples exploring such topics as evolution of altruism, competition between two inhomogeneous equations, a puts a new spin on the classical Lotka-Volterra predator-prey model. We also discuss the Fisher-Haldane-Wright equation, the Haldane principle for selection systems, and finally Fisher's fundamental theorem—the latter three topics we will return to again and again throughout the book.

In Chapter 6, we discuss models of frequency versus density-dependent population growth, and highlight some key differences between them. Chapter 7 dives into the discussion of inhomogeneous logistic and Gompertzian growth; we look at dynamics of distributions in inhomogeneous models and discuss various types of Darwinian (survival of the fittest) and non-Darwinian (survival of the common and survival of everybody) selection.

In Chapter 8, we introduce the Principle of minimal information gain, where we show that it can be derived from system dynamics rather than being postulated a priori. Chapter 9 discusses sub-exponential system dynamics and the Principle of minimal Tsallis information gain. The main result of this Chapter is that the Principle of minimal information gain is the underlying variational principle that governs replicator dynamics.

In Chapter 10, we discuss some philosophical issues on time perception and propose several hypotheses on how a model of inhomogeneous population extinction can be applied to time perception in a dying brain.

Chapter 11 explores several seemingly similar models of population growth—logistic, Gompertz and Verhulst, among others—and shows that intrinsic population composition may in fact be very different depending on which model described data best. We then apply this developed theory to cancer cell growth.

In Chapter 12, we show how the HKV-method can be applied to previously analyzed parametrically homogeneous systems to reveal the phenomenon of the expected value of the distributed parameter traveling through the phase-parameter portrait. This analysis can reveal new, complex and sometimes unexpected dynamical behaviors that help answer many interesting and important questions. We showcase several examples of this phenomenon, including the tragedy of the commons, natural selection between resource allocation strategies, and application of oncolytic virus therapy to a population of heterogeneous cancer cells.

Chapter 13 applies the HKV-method to game theory and looks at dynamics of selection of strategies within a single game. In Chapter 14, we look at selection between games, and then discuss how some of these insights can be applied to understanding the complex dynamics of cancer cells in tumors.

Finally, Chapter 15 (which can be read a stand-alone chapter), demonstrates how the HKV-method can be applied to selection systems with discrete time (maps).

Let us begin.

References☆

Bazykin A.D. Nonlinear Dynamics of Interacting Populations. World Scientific; 1998;vol. 11.

Bell G. Selection: The Mechanism of Evolution. Oxford University Press on Demand; 2008.

Berezovskaya F., Karev G., Snell T.W. Modeling the dynamics of natural rotifer populations: phase-parametric analysis. Ecol. Complex. 2005;2(4):395–409.

Brauer F., Castillo-Chavez C. Mathematical Models in Population Biology and Epidemiology. Springer; 2001;vol. 40.

Darwin, C. (1880). On the origin of species by means of natural selection: or the preservation of favoured races in the struggle for life. By Charles Darwin, … John Murray, Albemarle Street.

Dobzhansky T. Nothing in biology makes sense except in the light of evolution. Am. Biol. Teacher. 1973;35(3):125–129.

Hofbauer J., Sigmund K. Evolutionary Games and Population Dynamics. Cambridge University Press; 1998.

Johnson C. Introduction to Natural Selection. 1976.

Kuznetsov Y.A. Elements of Applied Bifurcation Theory. Springer Science & Business Media; 2013;vol. 112.

Nowak M.A., Sigmund K. Evolutionary dynamics of biological games. Science. 2004;303(5659):793–799.

Page S.E. Diversity and Complexity. Princeton University Press; 2010.

---

☆ To view the full reference list for the book, click here

Chapter 2

Inhomogeneous models of Malthusian type and the HKV method

Abstract

In this chapter, we present description and derivation of a new method for modeling population heterogeneity. We show why we refer to it as HKV, or hidden keystone variables, method and demonstrate how it allows incorporating a very high degree of heritable heterogeneity into dynamical systems. We then apply the method to several well-known growth models (namely, Malthusian, logistic, and Allee) and show how differently such populations can behave when population heterogeneity is taken into account (hint: we also see some non-Darwinian selection). The goal of this chapter is to provide initial exposure to the underlying mathematical theory, which will be deepened throughout the book, and to highlight that introducing population heterogeneity even into extremely well-studied models can reveal new, rich, and surprising dynamics.

Keywords

HKV method; Parametrically heterogeneous growth models; Non-Darwinian selection

2.1 Models of Malthusian type for inhomogeneous populations and a simplified version of the HKV method

Let us start from a simplified version of a model of inhomogeneous population growth and give two examples, which may clarify the main ideas of the developed approach. These examples provide motivation for focusing specifically on a wide class of models of evolving heterogeneous populations and corresponding replicator equations (RE). Furthermore, they provide a vivid illustration for how one can reduce a system of ODEs of initially large or even infinite dimensionality to low dimensionality through introducing auxiliary keystone variables (as was previously mentioned, the terminology was chosen to parallel the notion of keystone species in ecological systems, whose impact on overall system dynamics is disproportionally large relative to their abundance).

Consider an inhomogeneous population composed of individuals with different reproduction rates (Malthusian parameters) a; we refer to the set of all individuals with a given value of parameter a as an a-clone. Let l(t, a) be the size of a-clone at time moment t.

We assume that the growth rate of each clone depends on the total population size N(t). Dynamics of such a population can be described by the following model:

   (2.1)

where g(Nis the total population size. For example, if g(N) = const, then Eq. , then Eq. (2.1) describes an inhomogeneous logistic model, where each clone grows with an individual growth rate according to the logistic growth law up to the common carrying capacity C. (All of these and many other examples will be worked out in detail throughout the book after this mandatory introduction of the underlying theory.)

Now, let us derive expressions for the expected value of the distributed parameter a as it changes over time in the evolving population described by Eq. (2.1).

Denote the frequency of cell clone l(t, a; the probability density function (pdf) P(t, a) describes the distribution of parameter a throughout the population at each time moment tfor any integrable function f(a) to refer to the expected value of f(a) at time t. In a more general context described in (Hofbauer and Sigmund, 1998), it can be shown that the population size N(t) satisfies the equation

   (2.2)

and the pdf P(t, a) solves the replicator equation of the form

   (2.3)

where Et[a] is the mean value of parameter a.

Indeed, integrating Eq. (2.1) over a, we get

From this, we can calculate that

Now, let us assume that the initial pdf P(0, a) of the Malthusian parameter a is given and its moment-generating function (mgf)

   (2.4)

is known. This information could be obtained either from data when available or in the absence of data from general theoretical considerations.

In order to solve Eq. (2.1), let us define formally the keystone auxiliary variable q(t) as the solution to the Cauchy problem

   (2.5)

This equation cannot be solved at this moment, because the population size N(t) is unknown. However, clone densities and population size can be expressed with the help of the newly introduced in Eq. (2.5) keystone variable q(t):

   (2.6)

Therefore, by definition, total population size becomes

   (2.7)

Now the equation for the auxiliary variable q(t) can be written in a closed form:

   (2.8)

With this, we can completely solve the initial problem defined in Eq. (2.1) and corresponding replicator Eq. (2.3). Clone densities and population size are given by Eqs. (2.6), (2.7), respectively. More general results are summarized in the Reduction theorem, which will be proven in Chapter 4.

In what follows, we assume that the Cauchy problem defined in Eq. (2.8) has a unique solution q(t) in the interval 0 ≤ t < T, where 0 < T ≤ ∞.

Proposition 2.1

Current parameter distribution P(t, a) is determined by the formula

   (2.9)

The mgf of the current distribution P(t, a) is given by

   (2.10)

, one can see that Eq. (2.9) immediately follows from Eqs. (2.6), (2.7).

Next,

as desired.

It is well known that the moment-generating function can be used to compute all moments of a given probability distribution using the following formula:

where M(k)[0] denotes the kth derivative of the mgf M[λ] at the point λ = 0. The first and second moments, corresponding to expected value and variance of the distribution, will be of particular interest and value for our applications.

The mgf of the distribution P(t, a) is given by Eq. (2.10); hence

   (2.11)

The mean value Et[a] at each time moment can be then computed as

   (2.12)

The current variance can be computed as

   (2.13)

Proposition 2.2

The mean value Et[a] solves the equation

   (2.14)

Indeed,

Notice that Eq. (2.14) is a particular case of general Price’ equation (Robertson, 1968; Price, 1970, 1972); see also (Rice, 2006); this equation will be discussed in Chapter 5 within the framework of inhomogeneous population models.

We refer to Eq. (2.1) as a model of Malthusian type, and the variable q(t) can be viewed as internal time of the model, a kind of a measure of internal rate of evolution. With respect to the new time q, the model is composed of clones that grow as if they were independent from other clones. Indeed, making the change of variables dq → g(N)dt, we obtain from Eq. (2.6):

, then

The last equation is the standard Malthusian model. Hence, with respect to the internal time q(t), each clone grows according the Malthusian model with its own Malthusian parameter as it does not depend on any other clone or on the total population.

A more extensive discussion of the notion of internal time can be found in Chapter 10.

2.2 Dynamics of different initial distributions

Proposition 2.1 introduced in the previous section helps us trace the dynamics of the initial distribution of parameter a. Let us now consider some important examples.

Proposition 2.3

Let us assume that the initial distribution of the parameter a is

(i)normal with the mean a0 and variance σ0²; then, at any time moment, the current parameter distribution is also normal with the mean Et[a] = a0 + σ0²q(t) and with the same variance σ0²;

(ii)Poisson with the mean a0; then, at any time moment, the current parameter distribution is also Poisson with the mean Et[a] = a0eq(t).

Proof

Normal distribution

   (2.15)

has the mgf

   (2.16)

Then, according to Eq. (2.10),

It is the mgf of the normal distribution with the mean Et[a] = a0 + σ0²q(t) and variance σ0².

Poisson distribution

   (2.17)

has the mgf

   (2.18)

Then, according to Eq. (2.10),

It is the mgf of the Poisson distribution with the mean Et[a] = a0eq(t).□

As we will see soon, gamma distribution and its special case, the exponential distribution, are among the most important initial distributions for applications. Gamma distribution with coefficients k, s, η is defined by the formula

   (2.19)

a > η > 0, k > 0, and Γ(k) is the gamma function.

Gamma distribution has the mean

   (2.20)

and the mgf

   (2.21)

Proposition 2.4

Let us assume that the initial distribution of parameter a is Gamma distribution as given by Eq. (2.19). Let T∗ = inf {t : q(t) = s}. Then, parameter a is Γ-distributed at any time moment t < T⁎ with coefficients k, s − q(t), η such that

and

   (2.22)

Gamma distribution with η = 0, k = 1 becomes exponential distribution:

   (2.23)

.

It is useful to notice that the mgf of the exponential distribution can be written in the form

   (2.24)

Corollary

Let us assume that the initial distribution of parameter a is exponential as given by Eq. (2.23). Let T∗ = inf {t : q(t) = s}. Then, the distribution of parameter a is exponential at any time moment t < T⁎ with coefficient s − q(t) such that

   (2.25)

In applications, we will also consider versions of these distributions truncated on a bounded interval. In particular, exponential distribution truncated in the interval [0, b] has a form

   (2.26)

where a ∈ [0, bis a normalization constant.

The mgf of the probability distribution as given by Eq. (2.26) is

   (2.27)

Proposition 2.5

Let us assume that parameter a takes values in the interval [0, b] and its initial distribution is truncated exponential as given by Eq. (2.26). Then, at any instant, the parameter distribution is also truncated exponential in the same interval [0, b] with coefficient