When are population dynamics models useful?

When are population dynamics models useful? There seems to have been a lot of research about it, but how does it help? If I need data about how a population will evolve under what conditions, I need it because I need data for a decision (such as "can we kill 50% of population X without doing too much damage?"), right? But for that, the model needs to be aware of what causes what. And for that, I have to do experiments, right? Like "let's kill a significant amount of population X and see what happens in the next ten years". I really don't get it.

Population dynamics occupies a whole subset of mathematical biology. Perhaps the most pragmatic uses for modelling population dynamics come from the fields of epidemiology for modelling disease infection and transmission through a population (one such article), or ecology modelling things like forestation, fishing dynamics, predator-prey relationships (an example). Then there are more abstract uses, when you cannot measure a population to test a hypothesis because the labour is too intensive, or it's too costly, or no such surveillance mechanism exists. Theoretical models of population dynamics are built to gain an understanding of what the system as a whole may do under certain conditions. In this sense, population models are more of a theoretical exercise or a thought experiment.

Leonardo's already given you an excellent answer, but I thought I'd add my perspective. I'm a mathematical epidemiologist, so I'd at least like to believe these types of models are useful.

For me, there are a number of things population dynamics models are especially useful for:

• Highlighting data requirements. Yes, models need data, as you've mentioned. But they don't need all their data to come from one source, one study, or even one field. Models are also profoundly useful for showing where we don't have the data we need to fully understand a system. "To make a model where we understand A, we need the values for X, Y, and Z. X is well studied, but Y and Z aren't - though it turns out when we look over the entire parameter space for Z, nothing really changes in our answer. But guys? We could really use a study on Y."
• Eliminating guesswork. Models aren't perfect encapsulations of reality - there will always be some simplifying assumptions, etc. But it's better than "going with your gut" - especially for complex problems.
• Impossible studies. A ton of what mathematical epidemiology looks at it is areas where studies are either impossible, logistically difficult, or unethical. It would be very hard indeed to only be able to study pandemic response plans, or vaccination strategy only when we had an actual outbreak, or while a new vaccine is being rolled out.
• Highlighting potentially new directions. If you're considering an intervention, but no matter how effective you make it in your model it doesn't move the system much, it might not be worthwhile. Models can also highlight threshold effects - like the critical % of the population you'd need to vaccinate to achieve herd immunity.

Two previous answers listed many applications of population dynamics models. I want to add that they are also important for conservation of endangered species. For example classical stage-class model (Crouse et al 1987, free copy) indicate that the most effective way to protect sea turtles is reducing mortality of large juveniles.

Moreover, you don't have to perform such drastic experiments as killing 50% of a population to estimate your model parameters. Information about number of offspring, breeding success, natural mortality, etc can usually be gain without serious perturbation of wild populations. The number of individuals in every natural population fluctuates because of random reasons, so it is possible (but sometimes needs more field work) to calculate, (possibly nonlinear) regression between density of population and some demographic indicators and then extrapolate it to non-examined densities. For some small, short-living species it is sometimes possible to measure that correlations in laboratory or semi-wild conditions. For some long-living species, especially if they are sedentary, like trees, it may be better to compare specimens living in different distance from its neighborhoods. For poorly-known species it is possible to take missing information from related or similar species.

I'll throw one more application into the pot. Population dynamics also forms the foundations of population genetics, population ecology, and more recently plays an important role in frameworks such as evolutionary game theory and eco-evolutionary dynamics.

Here the models are also used as a type of theoretical exercise or thought experiment (as a previous answer suggests). In the development of evolutionary theory we simply cannot observe the process over the timescales we require to test the hypotheses we make. Thus, the development of population models allows us to explore 'possible worlds', as Robert May once put it, to see what kinds of adaptations or population structures we would expect to see, given the assumptions we put in.

We are also seeing an increasing number of population models and dynamical models used in conjunction with experiments on microbes in the field of experimental evolution. Here we can observe evolution in real time, and many of the assumptions about well-mixedness and large population sizes that are often made in modelling populations are actually fairly accurate.

Population dynamics models: Levins and the source-sink theory

For decades, ecologists and geneticists have tried to understand the population dynamics of different animal species, establishing models to understand the changes that occur in the number of individuals of a population, the composition of populations and the causes of such variations. Nature isn’t homogeneous or stable, so it’s not easy to set patterns to model these transformations.

Nature doesn’t work as a mathematical formula, in which populations inhabit a single area and can easily interact. For example, it’s unlikely for five males and five females of a species to be in the same area at the same time. In reality, landscapes are fragmented and these populations might live in different areas, isolated from each other. The fragmentation of a habitat hinders interactions among populations. By Mat Reding

In any case, population dynamics have tried to model the movements, changes and interactions that occur between groups of individuals, in order to understand the behaviour of different populations. In this article, we’re going to explain the foundations of some of these population models, which have been studied by different authors.

Lesson Plan

The following activity has been designed to teach basic principles of ecological population dynamics to high school and college level students. The activity utilizes online, interactive tools that allow students to manipulate characteristics of populations and analyze how these changes affect population sizes in real-time. This activity is best paired with a lecture on concepts such as (but not limited to) populations in ecology, carrying capacity, limited resources, or population growth. Students can complete the lesson with just a computer or tablet and internet access. Click here for the full lesson plan document.

Life Cycle of the Fruit Fly & the Population Modules

Fruit flies are used as a model organism in research laboratories and classrooms around the world. They are easily obtainable from suppliers, and many people have them in their own homes. They can be purchased cheaply, or they can be captured for free by placing a banana in an open jar for a day or two. Fruit flies are easily cared for, and there are no health or ethical concerns regarding their use in classrooms.

The life cycle of D. melanogaster is simple and predictable. Females lay up to 500 eggs on fermenting fruit. The male fertilizes the eggs, and the larvae emerge in 24–30 hours. The short, segmented, whitish-yellow larvae crawl around, gorging themselves on available food sources. They progress through three distinct phases within the larval stage, as denoted by changes in their mouth hooks (small, black hooks at their anterior it is not relevant to discern the different larval phases in the modules presented here). Toward the end of the larval stage, they begin climbing the walls of their enclosure, eventually becoming immobile as they transition into the pupal stage. At this time they turn dark brown and form a hard outer shell. They will remain in this state for about a week, after which time they emerge as adult fruit flies. An adult female fruit fly can begin to mate about two days after emerging.

Educational modules can be designed to explore a variety of biological processes using the life cycle of D. melanogaster as a model. Here, we use fruit flies to examine population dynamics. The first module explores population growth by comparing a population that begins with one female fruit fly and a population that begins with three females. Investigations focus on the total numbers of individuals within a population and the rates of growth across the different treatments. The second module examines the effect of food availability on population growth, and the third explores the effect of enclosure size on growth of the population. More specifically, this last module examines how the dynamics of similar populations differ if one population has significantly more space available than the other. Separately, each module highlights one limiting factor in population dynamics. Together, the modules illuminate relationships among different abiotic factors that influence population dynamics, with a specific focus on population growth. The three modules combined cover the three-dimensionality of NGSS (i.e., Scientific and Engineering Practices, Crosscutting Concepts, and Disciplinary Core Ideas) using authentic science investigations. Further, the modeling activities promote critical analysis and use of empirical data to construct working explanations in ways consistent with practices that Krajcik and Merritt (2012) recommended to promote student learning in science. NGSS connections for each module are presented in two tables at the end of this article.

Population Dynamics

Timothy D. Schowalter , in Insect Ecology (Fourth Edition) , 2016

Abstract

Populations of insects can fluctuate dramatically through time in response to changing environmental conditions. Disturbances are particularly important to population dynamics , triggering outbreaks of some species and locally exterminating others. Disturbances affect insect populations directly by killing intolerant individuals or indirectly by affecting abundance and suitability of resources or abundance and activity of predators and parasites. Population growth can be regulated (stabilized) to a large extent by density-dependent factors whose probability of effect on individuals increases as density increases and declines as density decreases. Primary density-dependent factors are intra- and interspecific competition and predation. Increasing competition for food (and other) resources as density increases leads to reduced natality and increased mortality and dispersal, eventually reducing density. Similarly, predation increases as prey density increases. Populations declining below their extinction threshold may be doomed to local extinction, whereas populations increasing above a release threshold continue to increase during an outbreak period. Development of population dynamics models has been useful for forecasting changes in insect abundance and effects on crop, range, and forest resources. General models include the logistic equation that describes a sigmoid curve that reaches an asymptote at carrying capacity. Chaos models have addressed the importance of initial conditions for subsequent changes in population size. More complex models incorporate specific population variables, including key factors and time lags. Despite limitations, models represent powerful tools for synthesizing information, identifying priorities for future research, and simulating population responses to future environmental conditions.

How Do Physical and Biological Factors Regulate Population Dynamics?

Patterns of population abundance are affected by a variety of biological and physical factors. For example, the abundance of a given species (for example, snails) might be controlled by the abundance of organisms that have a negative effect on the species of interest, such as competitors, predators, and diseases. Similarly, population abundance could be limited by the abundance of organisms that benefit the species of interest (for example, algae consumed by the snails).

In fact, some organisms require the presence of other species called symbionts with whom they live in direct contact. For example, corals use food molecules synthesized by symbiotic zooxanthellae (a type of algae), and zooxanthellae receive nutrients and protection from corals. However, not all populations are regulated by biological factors involving interactions with other species. Physical factors like water availability and temperature can control population abundance of some species.

Which type of factor (biological or physical) has a stronger effect on population dynamics? As one might suspect, the answer depends largely on the population that is studied. Some populations are regulated mostly by biological factors, others are controlled by physical factors, and most populations are affected by both biological and physical factors.

1.2 The laws of nature

The title of this chapter describes equation (1.1) as a “law” – what do I mean by law? Is it something that was out there, waiting to be discovered by humans, independent of our existence and thought? Or was it created by our thinking of it, consistent with reality but not of reality? Joe Rosen, formerly professor of physics at Tel Aviv University and the University of Central Arkansas, spent an entire volume thinking hard about these issues (Rosen 2010) . His categorization of reality and what we can know about it is useful and easy to follow, so I will use it here. He begins with the notion that there is an objective reality that exists independent of our existence. The primary reason for this observation is the simple fact that nature pushes back. Imagine a world where you can fly wouldn’t it be marvelous! If the world were not objective, but merely a construct of our imagination (a view of reality known as solipsism), then you could create this world, and fly. Unfortunately nature pushes back, and you will fall to the ground. So objective reality constrains what we can do.

The opposite of objective is subjective. Our inner thoughts and feelings are subjective, that is, they are known only to us as individuals. You might tell me what you are thinking or feeling, but I have no independent way of verifying that information. Beliefs about objective reality are similarly subjective, in that two people can have different beliefs about reality. However, it is possible for us to conduct reality checks on our beliefs about reality. If enough of us get together to check our beliefs, and over time, agree to a consensus belief that passes reality checks, then this is about as close to objective knowledge as we can get. Rosen calls this form of knowledge intersubjective it is different from subjective belief by virtue of its broader consensus amongst many people, and yet not fully objective by virtue of the fact that it was formed from our subjective perceptions of reality.

Intersubjective knowledge is socially constructed knowledge, but is not consistent with the post-modernist position that all reality is socially constructed. Our socially constructed, intersubjective beliefs are constrained by objective reality – not everything is possible. Even if a diehard post-modernist could convince a group of a 1000 people that she could fly unaided, she would not be able to do so.

In the field of wildlife management science, the goal is the production of “reliable knowledge” (Romesburg 1981) to use in making management decisions. It is not uncommon to see exhortations from leaders in the field to make science based decisions, presumably a call to use reliable, or intersubjective, knowledge to decide which course of action should be followed. Unfortunately, as we will see in many examples throughout this book, people do not make decisions like that. Our subjective beliefs about many things, from religion to justice, affect what we think should be done. Inevitably, the more people are affected by a decision or policy related to wildlife management, the more political (i.e. subjective) the decision or policy will become.

So a law of population dynamics is well tested intersubjective knowledge, or reliable knowledge, that we can use to make predictions about the consequences of management actions. As you will see below, a law will also have assumptions that must be met in order for it to apply.

4. Agent-based models and evolution

4.1. Population balance equations

Quantitative modeling plays an important role in bridging the relationship between single-cell and cell-population dynamics. The origins and consequences of cell-to-cell variability are often investigated analytically or computationally using single-cell models (e.g., [53, 91]). Such models typically ignore or idealize cell division and rarely capture the population-level effects of differential reproduction. An approach known as population balance modeling addresses this using partial differential equations [25, 92–97]. In this approach, cells are described as a continuous density flowing through a multi-dimensional state space that quantifies different physiological attributes (e.g., mass and chemical composition). Integral terms are used to account for birth and death processes along with a function describing the partitioning of cell contents at division. In the simplest one-dimensional case with no nutrient limitation or cell death, the population balance equation (PBE) can be formulated as follows

where F (x, t) is the number of cells, x is the cell mass, r (x) is the growth rate function of x, γ(x) is the division rate function that describes how the probability of a cell division varies with x, and is a partition function that describes the probability of a cell that divides with mass x to produce two sibling cells with masses x′ and x-x′. The first term on the left-hand side of Equation (33) is a transient term and the second is an advection term the right-hand side of the equation contains the source and sink terms.

The time-evolution is uniquely determined by the initial population distribution and all cell densities change according to the same deterministic rules. Correspondingly, information about individual cell trajectories is lost in this formalism. PBE models can be extended to include growth dependence on an external substrate and a diffusion term can be added to account for randomness in the evolution of cells in the state space. The modeling of discrete morphological stages or phases of the cell cycle can be represented by a set of coupled partial differential equations. Most PBEs cannot be solved analytically but can be discretized and integrated numerically. Unfortunately, when more details and higher dimensionality are incorporated, population balance models quickly become very difficult to formulate and solve computationally .

4.2. Cell growth and division

ODE models typically use first-order decay terms to account for dilution due to exponential growth . Similarly, discrete stochastic models approximate the effect of growth and partitioning of cell contents at division by increasing the degradation rates of all components. These methods tend to average away the dynamics resulting from growth and division, which can play an important role in cell dynamics (e.g., stochastic partitioning at cell division , asynchronously dividing cells , and asymmetric division [100–102]). In fact, it has been suggested by Huh et al.  that much of the cell-to-cell variability that has been attributed to gene expression “noise” (expression differences between genetically identical cells in the same environment) instead comes from random segregation at cell division, due to the difficulty in separating partitioning errors and gene expression noise profiles. Incorporating the details of cell division and gene expression into population models may help to resolve such controversies.

Computationally, cell growth and division can be modeled explicitly. Recent single-cell experiments show that the cell volume increases exponentially in bacteria [103–106], yeast [106–110], and mammalian cells [106, 111], and can be modeled using an exponential function [52, 112–114]

If cells grow at a rate that is proportional to the amount of protein they contain [115, 116], and if the protein concentration is constant, cells will grow exponentially in mass and volume . Modeling growth at the cellular level (see Section 4.4) can be important, as variability in single-cell growth rates may decrease the population growth rate . Such models for volume changes can be used to better capture intracellular protein dilution rates, since the dilution rate of any intracellular protein is given by d dilution = d dt In ( V ( t ) ) . This result can be derived by considering an intracellular protein at concentration c with copy number n at cell volume V, and where concentration of this protein drops only due to dilution resulting from cell growth. The rate of change of this protein is

If V (t) happens to depend on other cell types, then we can model dilution rates more generally using the above formula.

In addition to considering how the cell grows, one must also consider, based on the biology of the cell, how a cell regulates its size and how the cell volume and contents will be distributed at division. Cell size homeostasis may occur through different modes of regulation, including “sizers”, “adders”, and “timers”. Sizer regulation requires that a cell monitors its own size and cell division does not proceed until a minimal size has been reached . Adder regulation occurs when a cell adds a fixed size increment before division. Similarly, timer regulation describes the case when a cell grows for a fixed time duration before division. These three main modes of cell size regulation can be captured by the following equation 

Variability in size at division, or “sloppy cell-size control”, can be incorporated by treating cell division as a random process that takes place with a volume-dependent probability . Asymmetry in cell division can be modeled in either of these cases by setting the sum of the sibling cell volumes (V1 and V2) equal to the total volume of the dividing cell (V)

4.3. Constant-number Monte Carlo method

Models of non-interacting and non-dividing cells have been used extensively to study population variability arising from the process of gene expression. The population statistics from these models are often calculated from simulation of numerous independent realizations of individual cells [124, 125]. To incorporate cell division into the population’s dynamics, one might propose two approaches to investigate how populations evolve in time:

1. Simulate the time courses of single cells, randomly choosing one of the two newborn cells to follow when a cell divides. The result is lineages (or cell chains) containing a single individual per generation (e.g., ).

2. Simulate the time courses of single cells and continue to simulate all newborn cells produced. The result is a complete lineage tree (e.g., [122, 126]).

One problem with using the cell chain approach is that it does not account for the proliferative competition between cells in different physiological states, and so fails to provide the correct joint distribution of cell properties except in special cases. Hence, the second approach must be used when dealing with a model in which cell proliferation can vary with a number of intrinsic variables, such as age, metabolic state, and cell type. The problem here is that the size of the simulation ensemble rapidly grows to the point of intractability. This can be addressed using a Monte Carlo technique called the constant-number Monte Carlo (CNMC) method [127–129], originally developed to approximate the solution of population balance models (Section 4.1) of particulate processes.

The CNMC method is a statistically accurate method that keeps the total number of cells in an exponentially growing population fixed by randomly selecting cells at a user-specified time interval (which corresponds to experimental cell culture dilution times). This method is particularly suited to modeling well-mixed liquid cell cultures and has been applied to simulate heterogeneous cell populations [26, 52, 130]. The Gillespie algorithm [131, 132], which allows exact individual-based simulation of stochastic mass action kinetics, was combined with the CNMC method to accurately and efficiently simulate gene expression dynamics across growing and dividing cell populations .

4.4. Population dynamics algorithms

Population dynamics algorithms (PDAs) are computational frameworks that are important for studying cell population dynamics because they can account for a wide range of phenomena that cannot be investigated analytically or by simulating a model of a single cell. PDAs accomplish this by simulating a sufficiently large ensemble of individual cells, serving as a representative sample of the “true” population. This is advantageous because it puts the available methods for single-cell simulation at our disposal without the difficulty of integrating complicated and heterogeneous single-cell behavior into a broader mathematical framework. Notably, the use of individual-based models places virtually no constraints on the biologically relevant details that can be formulated and simulated . For instance, this allows biologically realistic features to be modeled with relative ease, such as cell growth and division effects (Section 4.2), that can be difficult to formulate and solve in an analytical framework.

Previous studies [123, 126, 134] paved the way for the individual-cell based PDAs. Though many of these previous approaches are extremely useful, they can be computationally prohibitive for simulating the dynamics of large, exponentially-growing cell populations. Motivated by this limitation, algorithms for the simulation of intracellular content and cell growth and division were developed, ranging from deterministic  and stochastic Langevin approaches , along with methods to determine the timing of cell divisions and partitioning of cell contents that agree with population balance formulism (Section 4.1), to parallel stochastic approaches that described growing and dividing cells [52, 114]. One method, which we refer to here as the “asynchronous PDA”, combines the Gillespie stochastic simulation algorithm [131, 132] with a CNMC method (Section 4.3) to simulate population dynamics in a computationally efficient manner (Fig. 3A). Here, “asynchronous” refers to cells being simulated independently from each other and the global population state being determined (and statistics calculated) at discrete and usually equally spaced synchronization barriers to permit parallelization of the algorithm. In the asynchronous PDA, the population is restored at the synchronization barriers to a predefined size using the CNMC method. An accelerated version of the asynchronous PDA, which is applicable when steady-state and symmetric cell division can be assumed, was subsequently developed . The accelerated PDA is well suited for expediting simulations by performing coarse-grained explorations of parameter space, to be subsequently investigated in more detail using the asynchronous PDA.

Population dynamics algorithms. (A) Flow diagram for the asynchronous population algorithm, where all cells are simulated independently of one another and synchronized only when the simulation time for each cell (ti) is equal to or exceeds the user specified sampling time (tsample). The population size is restored using a constant-number Monte Carlo (CNMC) method to a prespecified fixed size each time the simulation time (t) is greater or equal to the population restore time (trestore). Xi is the system of equations/reactions and Fi the fitness that correspond to cell i, respectively. (B) Schematic illustrating the concept of a general population simulation framework. Which reaction occurs next (Ri) and the time at which it will occur (ti) in each cell is determined stochastically. This approach allows for intracellular communication (represented by purple triangles and arrows) and resource consumption (represented by blue squares and arrows) in “real-time” [as opposed to only at each tsample in (A)]. Here, the next reaction will occur for cell 1 (i1 = R1) at t1 = 1.12, when it will uptake a signaling molecule exported from cell 3 at an earlier time. Fortran code for (A) is available in the Appendix B of  and at: https://github.com/dacharle/PDA_Fortran, and an object-oriented C++ prototype at: https://github.com/alanyuchenhou/gene-expression (Color online). Another more general framework, inspired by the concept of “reaction channels” in Gillespie’s algorithm [131, 132], links simulation channels through scheduling dependency graphs (introduced by Gibson et al.  to improve the performance of the Gillespie algorithm) to handle the scheduling and execution of state-updating events on individual cells . This approach is analogous to Gillespie’s algorithm, where the propensities of different reaction channels are used to determine when the next reaction will occur (scheduling) and how the numbers of molecules are changed when it occurs (execution) (Fig. 3B). Simulations are performed using an asynchronous method, which is ideal and parallelizable for non-interacting cells, and a synchronous method, enabling the incorporation of cell-to-cell communication (which is not practical in the PDAs due to the parallel nature of the algorithms). For a discussion of synchronous versus asynchronous models, see ref. . Once the population size limit is reached, the CNMC method (Section 4.3) is introduced (as in the asynchronous PDA) to keep a fixed sample population size with the appropriate composition.

In summary, the dynamics of heterogeneous cell populations can be highly complex and are difficult to investigate analytically. The frameworks presented in this section address this by enabling efficient individual-based population-level simulations without the need to formulate or solve complex mathematical equations. They are designed specifically to ease the incorporation of user-designed biological features and to facilitate the transition towards population-level modeling in quantitative biology.

4.5. In silico evolution

Many of the mathematical models and computer algorithms discussed so far in this review article can be modified to model evolution. This is important, for example, because it allows us to perform long-term in silico evolution experiments in scenarios that may be difficult or costly to investigate in the laboratory. In Section 4.5.2 we present a simple computational model of evolution more complex computational frameworks to model the evolution of a cell population are discussed in Section 4.5.3.

4.5.1. Fitness

Darwin’s theory of evolution by natural selection is built on the idea that some genotypes have higher fitness than others. However, what exactly mean by the term “fitness” is not always clear and term has been used to mean subtly different things . In fact, even the unit of selection, whether it be the gene or individual [137, 138] or group (recently rebranded as multilevel selection theory) [139, 140], is still debated. A related concept is inclusive fitness, which describes the total effect an individual has on proliferating its genes by producing off spring and by providing aid that enables relatives to reproduce . An evolutionary strategy for increasing inclusive fitness, even at a cost to the individual’s own survival and reproduction, is known as kin selection. According to Hamilton’s rule, kin selection causes genes to increase in frequency when the genetic relatedness of a recipient to an actor multiplied by the benefit to the recipient is greater than the reproductive cost to the actor. Fitness landscapes have long been used to illustrate the effect of genetic factors on fitness , and more recently, nongenetic factors as well [8, 13].

The exponential growth rate (Section 2.1) is one common measure of fitness in microbiology and experimental evolution studies. At the population level, this is done by measuring the number of cells or the optical density of the cell culture (e.g., using a cell counter or a spectrophotometer, respectively) and then obtaining the growth rate by fitting the data to an exponential function [Equation (2)]. Population fitness is also commonly measured by direct head-to-head competition assays . This is often done experimentally by determining the relative fitness of each competitor with respect to a reference strain. For example, non-fluorescing evolved cells can be competed (and distinguished) against an ancestral strain that expresses the green fluorescent protein . The fitness (W) is then calculated by

Cell growth rates within a clonal population can vary depending on the environmental context in which it evolved. While a constant environment selects for low variance in growth rate, a fluctuating environment can select for high variance if the growth rate correlates with survival under stress . Thus, the growth-rate distribution is an important evolutionary parameter that can be captured using single-cell experimental measurement techniques  and modeled using population simulation algorithms (Section 4.4).

Population fitness is distinct from the cellular fitness of its constituent members, though the former can be obtained from the latter. For instance, we can define the population or “macroscopic” fitness W of an isogenic population under stress as 

4.5.2. Ordinary differential equation evolution model

The model presented in this section describes how the number of cells with wild-type and mutant genotypes varies over time based on their fitness . Specifically, this model describes population dynamics by a system of ordinary differential equations and assumes a constant population size and mutation rate. Here, wild-type and mutant cells are characterized by a single fitness parameter. This ODE evolution model (a corresponding but more detailed evolutionary simulation framework is discussed in Section 4.5.3) has three free parameters: rate of beneficial mutations μ, input probabilities P(G) and P(T) of a given mutation being type G (genomic) or, K (knockout) or T (tweaking) P(K) is determined via P (K) = 1 - P (G) - P (T). Note that while the probability of P(K) could be 1, its rate μP (K) is ⪡1 per genome per generation.

The approach taken in this model is similar to that presented in Section 4.1, where the “gain” and “loss” rates of each genotype are used to develop a system of ODEs that describes the population size of each genotype i over time. Assuming that the number of beneficial mutations arising per unit time is proportional to the number of cell divisions (i.e., mutations arise strictly due to DNA replication errors), it can be described by

If we consider only the influx of new genotypes that survive drift (i.e., assuming all other mutations go extinct rapidly), then the effective influx of genotypes Mi that carry a potentially beneficial mutation of type i equals

This model can now be applied to specific cell types by defining the ancestral/mutant types (M0/Mi) and the associated fitnesses (f0/fi). This modeling approach is more general than the more detailed computational approach that is briefly discussed in the following section and facilitates large-scale parameter scans.

The ODE evolution model predicted how fast experimental wild-type genotype disappears from the population, as well as the mutation type (K, T, G) that predominantly replaces the wild type in each condition . Interestingly, the ancestor genotype disappeared fastest in conditions with steep monotone cell fitness landscapes (see  for a review of fitness landscapes) and remained in the population longer in peaked cell fitness landscapes each experimental condition favored different fractions of mutations types as long as they were available.

4.5.3. Evolutionary simulation frameworks

More detailed computational evolutionary simulation frameworks to model molecular evolution have been developed that explicitly account for experimental details (such as phenotypic switching and resuspensions). One such framework from the same study  as the ODE evolution model described in Section 4.5.2 enables the prediction of how experimental evolution will affect evolutionary dynamics (Python code is available in the supplemental materials of Gonzalez et al. ). A more detailed computational framework was required because the simpler ODE evolution model could not predict the number of distinct mutant alleles in the evolving population and lacked important experimental details (e.g., periodic resuspensions and phenotype switching between ON and OFF states with experimentally determined switching rates). The evolutionary simulation framework predicted the number of distinct mutant alleles, in addition to the characteristics predicted by the simpler mathematical model. Importantly, the modeling in Gonzalez et al.  was crucial for understanding the evolutionary dynamics of mutants arising, establishing, and competing, as well as the number of alleles in the cell population.

Another computational framework that can be used to model the evolution of cell populations is the asynchronous PDA (Section 4.4), which was recently modified to incorporate evolution . This was done by modeling genetic mutation as a change in the reaction rate parameters of a mutant cell probabilistically each time a cell divides. As cell fitness (cell cycle time) is coupled to gene expression level and selection pressures, this allows for selection of the most fit genotype. Importantly, this approach provides a framework for studying how nongenetic variability in gene expression can affect evolution and predicted that the level of evolved cell-to-cell variability in the population depends on the associated fitness costs and benefits of gene expression in a specific environment. An exact algorithm for fast stochastic simulations of evolutionary dynamics was developed by Mather et al., , which provides a significant speedup when the population size is large and mutation rates are much smaller than the birth and death rates.

INTRODUCTION

The spatial distribution of individuals in a plant population is mainly determined by seed dispersal patterns and subsequent establishment success. Most wind-dispersed seeds land near their source and very few travel over large distances (Bullock and Clarke, 2000 Nathan and Muller-Landau, 2000 Tackenberg, 2003). Seedling establishment in orchids can be facilitated by conspecific plants because the chance of forming associations with beneficial mycorrhizal fungi is higher close to established plants (Diez, 2007). However, with increasing plant density, fecundity and survival in a population decrease as a result of increasing competition (e.g. Watkinson and Harper, 1978 Mithen et al., 1984), or as a result of increased attraction of herbivores and seed predators (Feeny, 1976). Dispersal, in turn, is not only influenced by species-specific seed characteristics, but also by the distribution of adult plants and local demographic processes, primarily fecundity (Clark and Ji, 1995 Levin et al., 2003).

In a metapopulation, seed dispersal is necessary to colonize unoccupied habitats, to re-colonize extinct patches, and to transfer seeds from high to low competition patches (Husband and Barrett, 1996). However, high dispersal distances reduce local fecundity and the probability that seeds will reach a suitable habitat. Thus, there is a conflict between seed survival and colonization (Levin et al., 2003), the existence of which was demonstrated for a hypothetical metapopulation (Johst et al., 2002) and for Pinus halepensis (Bohrer et al., 2005). Whether long-distance dispersal is beneficial depends on local population dynamics. For instance, long-range dispersal has a positive effect on metapopulation persistence unless the number of potential dispersers is low due to small population growth rates (Johst et al., 2002). The interactions between fecundity and mean dispersal distance, and the consequences for population maintenance, are little studied.

The effect of spatial distribution in populations is now an essential component of metapopulation models, and has become firmly established in population biology (Hanski and Simberloff, 1997). Classical metapopulation theory focuses on the extinction and re-colonization of local patches regarding neither size and spatial arrangement of the patches nor local population dynamics (Levins, 1969). More recently, metapopulation theory has been considerably extended to include spatially explicit and realistic approaches including declining non-equilibrium, habitat-tracking, and mainland-island metapopulations (Thomas, 1994 Hanski, 1997 Harrison and Taylor, 1997). Most metapopulation models, however, concentrate on regional dynamics and ignore details on the scale of local populations (but see Bohrer et al., 2005 Volis et al., 2005 Mildén et al., 2006).

Previous studies on the population dynamics of vascular epiphytes found survival to be the most important parameter determining population growth rates (Hernández-Apolinar, 1992 Zotz et al., 2005 Zotz and Schmidt, 2006 Winkler et al., 2007), but none of these studies accounted for the spatial structure of epiphyte populations. In the present study, it is demonstrated how the incorporation of spatial structure alters model predictions of population size and the importance of demographic processes for population growth rates of three epiphytic orchids. Spatially realistic matrix metapopulation models are used, where population dynamics at the scale of the local population (individuals growing on one host tree) are based on detailed stage-structured observations of transition probabilities and epiphyte populations on different trees are connected by a dispersal function. The question is asked whether demographic processes differ significantly between trees, and if so, what the consequences for metapopulation structure and persistence are. Furthermore, metapopulations where increased dispersal results in reduced local fecundity are compared with hypothetical models where local fecundity is not related to dispersal.

Point 1. Model Population Dynamics

The Wilson-Cowan Model

The Wilson-Cowan model represents the dynamics between the spatially confined excitatory and inhibitory subpopulations. Here, (E(t)) and (I(t)) represent the instantaneous discharge rate of the excitatory and inhibitory population at time (t) , respectively. We define the variables characterizing the dynamics of a spatially localized neural population as the following.

Now, we have the equations (E(t)) and (I(t)) . We assume that the value of these functions at time ((t+ au)) will be equal to the proportion of cells, receiving at least threshold excitement at time (t) , that isn’t sensitive (not refractory).

Neuron Assumptions

Assumptions. We are concerned with the behavior of the subpopulations rather than individual cells. Under the following assumptions, we may ignore random spatial interactions:

1. Assume that the cells containing the two populations are very close in space
2. Interconnections are arbitrary yet dense enough so that it is likely there’s at least one path connecting any two cells in the population
3. Neural processes depend on the interaction of excitatory and inhibitory cells

Support material.

1. Proportion of Sensitive Cells in each sub population

First, the study obtained the independent expression of the proportion of sensitive cells and the proportion of cells that received at least a threshold stimulation. Relevant biologic assumptions assess when a neuron fires. Consequently, the most simple model must fulfill two conditions:

1. The neuron should not be in its “refractory period,” that is, it can’t fire again immediately after having fired.
2. Neurons need to receive enough input in a short time.

If the total refractory period has a duration of (r) at time (t) , then the proportion of excitatory cells that are refractory will be given by the following:

Therefore, if (r) is the length of the refractory period, the proportion of excitatory cells that are sensitive is the fraction of neurons that satisfy condition one at time (t) :

Similar expressions are obtained for the inhibitory subpopulation.

2. Subpopulation Response Functions

If all cells receive same number of excitatory and inhibitory afferents, then the subpopulation response function (S(x)) will take the form:

where (D( heta)) characterizes the distribution of individual neural thresholds.

In other words, we assume that all cells in a subpopulation have the same threshold ( heta) and that there is a combination of the number of afferent synapses per cell. Suppose (C(w)) is the synapse distribution function, and (x(t)) is the average excitation of each synapse. In that case, we expect all cells with at least ( heta) (x(t)) synapses to receive sufficient excitation.

If all cells within a subpopulation have the same threshold ( heta) , then the distribution of the number of afferent synapses per cell characterized by (C(w)) . Therefore, the subpopulation response function takes the following form:

The study assumes that the resting threshold ( heta) is the same for all cells in the population. Hence, it holds that the sigmoid response function is related to the distribution of synapses.

Sigmoid Response Function

Population Growth Models. The sigmoid, logistic growth model is commonly used in biology to model the growth of a population under density dependence:  Figure 1: Logistic Growth Model. Plot of the flow field, horizontal lines at any equilibrium points, several trajectories for the case (eta=1) and (K=4) .

To represent neurons’ nonlinear behavior, we use a sigmoid function (S(x)) , taken as characteristic of any subpopulation response function, dependent on two parameters (a) and ( heta) :

The response function is sigmoidal if (S(x)) monotonically increases on the interval (- (infty) , + (infty) ), approaches the asymptotic states zero and one as (x) approaches – (infty) and + (infty) respectively, and has one inflection point. Figure 2: Sigmoid Subpopulation Response Function. The particular function shown here is the logistic curve: (S(x)=frac<1><1+e^<-a(x- heta)>>,) such that ( heta=5, a=1) . (X) is average level of excitation in threshold units.

Sufficient excitation. To see if condition 2 is fulfilled, we need the total input to the subpopulation to be a weighed contribution:

We call (S_e) the response function, also called the input-frequency characteristic of the excitatory neurons, and correspondingly (S_i) for the inhibitory ones. A spike can still help eliciting a new spike in a downstream neuron even a few milliseconds later. The response functions give the expected proportion of those receiving at least threshold stimulation as a function of the population’s average excitation levels. The sigmoid, nonlinear functions (S_e) , and (S_i) take the form: Figure 3: Plot of the Sigmoid functions (S_e) and (S_i) for the excitatory and inhibitory subpopulations, respectively. Parameters: (a_e = 1.3, heta_e=4, a_i=2, ext < and > heta_i=3.7)

3. Average Level of Excitation

The inhibitory potential is when hyperpolarization brings about a net negative charge, so the potential is further away from zero. The excitation potential occurs during the depolarization process, so the potential is closer to the excitation threshold. Here is an expression for the average level of excitation generated in the cell at time t:

[ egin int^t_ <-infty>alpha (t-t^prime) [c_1E (t^prime) - c_2I (t^prime) + P(t^prime)] dt^prime end ag <7>]

Here, we define the parameters (c_1, c_2>0) as the connection coefficients, representing the average number of excitatory and inhibitory synapses per cell. We denote (P(t)) and (Q(t)) as the external input of excitatory and inhibitory subgroups, respectively. The inhibitory subgroup will use similar expressions but with different coefficients and different external inputs.

The different coefficients reflect the difference in axon and dendritic geometry between excitatory and inhibitory cell types. In contrast, variations in external inputs assume the existence of cell types specific to the population. We can illustrate the diversity in axon and dendritic geometry between excitatory and inhibitory cell types by drawing several neurons, as shown below.    Figure 4: Axonal and Dendritic geometry. Plot of Several Neurons to illustrate the diversity in axonal and dendritic geometry between excitatory and inhibitory cell types.

Overall, the dendrite branch structure is an essential feature in coupling synaptic inputs and managing action potentials in neurons. Each neuron has a unique branch density and pattern: this unique morphology correlates to the neuron’s function.

4. Equations for the activities E(t) and I(t)

Suppose the probability that a cell is sensitive is independent of the probability that a cell is excited above the threshold. In that case, we can define the excitatory subpopulation as the following equation:

We designate this correlation between excitation and sensitivity by the following expression (gamma left[E(t^prime) dt^prime,mathcal_e(x)> ight],) so that the previous expression becomes:

[ small egin &left[ 1-int^t_E(t^prime) dt^prime ight] &qquad qquad > mathcal_e(x) left<1-gamma left[E(t^prime) dt^prime, mathcal_e(x)> ight] ight > delta t end ag <9>]

Due to fluctuations inherent in average excitation & cell thresholds, the parameter (gamma) is taken to be zero: (gamma = 0)

5. Dynamics of a Localized Population of Neurons

It follows that the equations governing the dynamics of a localized population of neurons take the form of the following expressions:

[ egin mathbf(t+ au) =& E(t^prime) mathrmt^prime ight ]> &quad > mathcal_e left < int^t_<-infty>alpha(t-t^prime) [c_1 E(t^prime) - c_2 I(t^prime) + P(t^prime)] mathrmt^prime ight> , ag <10>end ]

for the excitatory and inhibitory subpopulations. The parameters ( au) and ( au^) represent the response delays, after which cells at time (t) will be firing.

6. Time Coarse Graining

We are interested in the coarse-grained transitory behavior of the neural activity. Coarse-grained modeling aims to simulate the behavior of complex systems using simplified representations. In the above equation, we can get rid of the integral and multiply the stimulus by a constant describing the length of time influence. Consider the following application of time course-graining of equation E(t):

Now, assuming that (alpha(t-t^)) decays exponentially and excludes inhibitory interactions, we extend the equation to the lowest term in ( au) :

Now, we perform a Taylor series expansion about ( au=0) , which is the retention of the linear term:

We can say that the activity at time (d + dt) depends on the simultaneous satisfaction of conditions (1) and (2):

[ E(t + dt) = left(1- rEleft(t ight) ight), S_e , left(kc_1 Eleft(t ight) - kc_2 Ileft(t ight) + kPleft(t ight) ight) ] ,

For physiological significant values of (alpha) and (r) , the course grained equations are valid.

Transition and Connection

After the above reductions and assumptions, turning the original system in differential form and suitably rescaling (S_e) and (S_i) , we reach a coupled, nonlinear, differential equation for the excitatory and inhibitory populations’ firing rates.

[ au_e frac

= -E + (k_e - r_e E) , S_e(c_1 E - c_2 I + P)]

[ au_i frac

= -I + (k_i - r_i I) , S_i(c_3 E - c_4 I + Q),]

where we define ( au_e) and ( au_i) as the time constants, (k_e) and (k_i) as the non-dimensional constants, (r_e) and (r_i) as the constants describing the length of the refractory periods, (S_e) and (S_i) as sigmoid functions representing the nonlinearity of the interactions, (c_<1,2,3,4>) as the parameters representing the strength of the excitatory to inhibitory interactions, and (P) and (Q) as the external inputs to the excitatory and inhibitory populations respectively.

The System’s Resting State. The point ((E=0, I=0)) is the system’s resting state, which we require to be a stable fixed point. The resting state’s mathematical result is that (E=0, I=0) must be a steady-state solution for (P(t)=Q(t)=0) , i.e., in the vacancy of external inputs. This can be accomplished by transforming (S_e) and (S_i) so that (S_e(0)=0) and (S_i(0)=0) . Thus, we subtract (S(0)) from the original function. We have that the maximum values of (S_e) and (S_i) are less than 1, which we denote as (k_e) and (k_i) . Furthermore, the resting state must be stable to be of physiological significance.

[ egin au_efrac

&= -E + (k_e - r_e E)S_e(c_1 E - c_2 I + P), ag <12> au_ifrac
&= -I + (k_i - r_i I)S_i(c_3 E - c_4 I + Q) ag <13>end ]

We can write the equations for the nullclines corresponding to (dE/dt=0) and (dI/dt=0) . Since (S_e) and (S_i) are both sigmoidal, the functions have unique inverses. Hence, we can equate the nullclines to the following expressions:

[ egin c_2 mathbf &= c_1 mathbf + S_e^<-1>left( frac ight) +P ,quad frac

=0 ag <14> c_3 mathbf &= c_4 mathbf + S_i^<-1>left( frac ight) - Q ,quad frac
=0 ag <15>end ]

Integrating genomics in population models to forecast translocation success

Department of Fish and Wildlife Sciences, University of Idaho, Moscow, ID, U.S.A.

Address correspondence to T. Seaborn, email [email protected]

Institute for Bioinformatics and Evolutionary Studies (IBEST), University of Idaho, Moscow, ID, U.S.A.

Biological Sciences, Boise State University, Boise, ID, U.S.A.

Department of Biological Sciences, Idaho State University, Pocatello, ID, U.S.A.

Department of Fish and Wildlife Sciences, University of Idaho, Moscow, ID, U.S.A.

Biological Sciences, Boise State University, Boise, ID, U.S.A.

Biological Sciences, Boise State University, Boise, ID, U.S.A.

Department of Fish and Wildlife Sciences, University of Idaho, Moscow, ID, U.S.A.

Address correspondence to T. Seaborn, email [email protected]

Institute for Bioinformatics and Evolutionary Studies (IBEST), University of Idaho, Moscow, ID, U.S.A.

Biological Sciences, Boise State University, Boise, ID, U.S.A.

Department of Biological Sciences, Idaho State University, Pocatello, ID, U.S.A.

Department of Fish and Wildlife Sciences, University of Idaho, Moscow, ID, U.S.A.

Biological Sciences, Boise State University, Boise, ID, U.S.A.

Biological Sciences, Boise State University, Boise, ID, U.S.A.

Author contributions: all authors contributed to reviewing papers, writing the review, creating the figures, and editing the manuscript.

Abstract

Whole-genome sequencing is revolutionizing our understanding of organismal biology, including adaptations likely to influence demographic performance in different environments. Excitement over the potential of genomics to inform population dynamics has prompted multiple conservation applications, including genomics-based decision-making for translocation efforts. Despite interest in applying genomics to improve translocations, there is a critical research gap: we lack an understanding of how genomic differences translate into population dynamics in the real world. We review how genomics and genetics data could be used to inform organismal performance, including examples of how adaptive and neutral loci have been quantified in a translocation context, and future applications. Next, we discuss three main drivers of population dynamics: demographic structure, spatial barriers to movement, and introgression, and their consequences for translocations informed by genomic data. Finally, we provide a practical guide to different types of models, including size-structured and spatial models, that could be modified to include genomics data. We then propose a framework to improve translocation success by repeatedly developing, selecting, and validating forecasting models. By integrating lab-based and field-collected data with model-driven research, our iterative framework could address long-standing challenges in restoration ecology, such as when selecting locally adapted genotypes will aid translocation of plants and animals. During this time of mass disruption, be advised that we appreciate there will be a slower pace for all. Restoration Ecology understands that reviews and decisions may be delayed responses from authors may be delayed. There are no consequences for delays. We ask all to be patient. The EIC and Managing Editor work remotely as is (in different countries) so we already work from ‘home’.

We are attempting to add this message to our communications (not as easy because the Editors don’t have total editing rights) and reduce the normal reminder emails to reflect this uncertain time. If you receive our normal email correspondence reminding you of deadlines, we are waiving these and asking only that you let us know, if possible, of delays exceeding a month.