SexOnR–04_Selection. Four aspects of selection and an example of simulating differential mortality
The creation of the first model in R is described in detail. For now without sex.
4.2 Simulating selection4.2.1 Selection, its four aspects and four possible stages of its simulation
What is natural selection? Differentiated, i.e., dependent on the properties of the units of selection, their preservation and reproduction. In the most straightforward case (which is far from being the only possible one), for organisms, selection is their differentiated, trait‑dependent survival and reproduction.
Selection should not be regarded as a purely biological phenomenon. The significance of selection was first consciously recognized in biology. This awareness is most often associated with the name Charles Darwin (Charles Darwin, 1809–1882), although a substantial contribution to this awareness was also made by other brilliant British scientists Thomas Malthus (Thomas Malthus, 1766–1834), Alfred Wallace (Alfred Wallace, 1823–1913) and Herbert Spencer (Herbert Spencer, 1820–1903). Nevertheless, selection in biology is a particular case of a much broader phenomenon. Contemporary American philosopher Daniel Dennett (Daniel Dennett; b. 1942) considers the idea of selection the most important idea of humanity. Thanks to the work of the great Austrian physicist Ludwig Boltzmann (Ludwig Boltzmann, 1844–1906) it became clear that in nature processes occur spontaneously that lead to more stable states. Stability is the ability of a given state to persist (and possibly to spread). According to this approach, the architect of the universe is selection sensu lato (in the broad sense): more stable processes, i.e., those that with higher probability preserve, restore or spread their state, displace less stable processes. The process of survival and preservation of better‑adapted organisms (i.e., those whose traits match their environment better and increase the probability of their survival and reproduction) is more stable than that of less adapted organisms. We live in a universe that is shifting toward more stable processes, among which is the existence of increasingly well‑adapted organisms.
In our models we will simulate selection mainly at the level of organisms, occasionally at the level of populations (group selection), possibly at the level of genes. One way or another, selection phenomena will be observed either in differentiated survival or in differentiated reproduction. Different reproductive probabilities constitute a separate form of selection, sexual selection, which was discovered by Charles Darwin. Note: in simulation models (i.e., models that describe processes in the studied systems “step by step”), survival and reproduction will most likely be simulated separately, as sequential processes. In reality this is indeed the case: an organism must first reach maturity before it can reproduce. Thus, survival and reproduction become distinct stages of the work cycle of a typical simulation model in our course. Moreover, in the typical case we will simulate reproduction as a two‑step process. To understand why we make this choice, one must understand what survival depends on.
In reality living organisms interact with a huge number of environmental factors. On the one hand there is the abiotic environment, the biotope. One must obtain resources from this environment (water, air, light, various minerals) and protect oneself from its adverse conditions (cold, heat, drought, toxic substances, etc.). There are also myriad interactions with the biotic environment, the biocenosis. There exist several classifications of relationships with other organisms, but in this context the most important distinction is between individuals of other species and conspecifics, members of one’s own population.
In the textbook by D. Shabanov and M. Kravchenko “Ecology: the biology of interaction” it is proposed to divide the environment as a whole into sub‑environments: the biotopic sub‑environment (the set of important abiotic factors related to the place of residence), the biocenotic sub‑environment (relationships with organisms of other species, predators, prey, parasites, symbionts, etc.) and the population sub‑environment (the set of relationships with other individuals in the population). In humans additional sub‑environments are added, but they are not relevant for our discussion. What is important for us is to distinguish two types of mortality: one not directly linked to population density (or size), and one linked to intra‑population competition, which leads to a reduction in density or size.
What happens when the “Titanic” sinks? That ship carried a large number of people differing in sex, age and health status. These differences meant that different people had different probabilities of dying from different causes. Mortality on that ship was not high, but some could die in accidents, some could be poisoned by alcohol or poor food, some could die of disease. Then the whole passenger body of the “Titanic” (which can be regarded as a population defined by a formal characteristic) underwent a sharp reduction in number, essentially down to the number of seats in the lifeboats. Among the passengers with the cheapest tickets, mostly healthy men survived because they were able to take places in the boats; among first‑class passengers (different classes were rescued from different decks) mostly children, women and the elderly were let ahead of men. In any case, the pattern of mortality that was linked to intra‑population competition (during boarding the lifeboats) was completely different from that which was not linked to it.
Adaptations for survival in the external environment and adaptations for survival under intra‑population competition can be different. In reality sometimes death and survival occur randomly, i.e., the fate of each organism does not depend on its traits. Simulation models can include such phenomena, but they are not examples of selection. Success or failure in reproduction is another “degree of freedom”. Here we can see that reproductive success depends, at minimum, on two variables. The first is the opportunity to participate in reproduction (e.g., to form a pair). The second is the expected number of offspring in that pair. These variables can be relatively independent.
In the Siverskyi Donets basin of the Kharkiv region, hemiclonal population systems of the hybridogenic complex of green frogs are characteristic, in which the breeding assemblage includes individuals of the species Pelophylax ridibundus (Pallas, 1771) as well as interspecific hybrids between Pelophylax ridibundus and Pelophylax lessonae (Camerano, 1882). These hybrids are known under a separate (species‑like) name, Pelophylax esculentus (Linnaeus, 1758). Males of P. esculentus behave on the breeding ground much more aggressively than males of P. ridibundus; hybrid males are able to repel females of the parental species. Thus, P. esculentus has a higher chance of forming a pair. However, the mechanism of gamete formation in such hybrids is considerably complicated, because they typically produce not hybrid gametes but gametes belonging to one or the other parental species. In many individuals this gamete formation is disrupted, resulting in a large proportion of eggs laid with such a male remaining unfertilized. The number of offspring of a P. esculentus individual that formed a pair is often much lower than that of a P. ridibundus male that succeeded in obtaining a female. The fecundity of P. ridibundus is higher.
Thus, in simulation models that study processes in populations, it may be appropriate to simulate four aspects of selection:
— differentiated survival and death; survival probability should be proportional to viability, and death probability to mortality;
— competitive reduction of numbers; the probability of persistence under such reduction is proportional to competitiveness, and the probability of loss to competitive vulnerability;
— differentiated participation in reproduction; for individuals that must find a mate, the probability of successful pair formation can be measured by attractiveness;
— differentiated fecundity; different organisms can be described by different fecundities.
In this section we will examine various ways to simulate differentiated mortality, while the simulation of the other selection aspects will be considered later.
Population <- 1:10
Population
## [1] 1 2 3 4 5 6 7 8 9 10
Mortality <- sample(0:10, 10, replace=T)/10
Mortality
## [1] 0.1 0.2 1.0 0.3 0.7 0.8 0.8 0.9 0.2 0.1It is important for further discussion that the Population vector lists individual codes, the Mortality vector lists the probabilities of death for individuals at the same position in the Population vector, and these vectors have the same length. If all individuals in the Population vector belonged to the same form and were subject to selection of the same strength, the Mortality vector would be filled with identical data. We are considering a more complex case. For the modeling method discussed, there is no difference between these cases. In section 4.1.2, we described the first modeling path, simulating individual "fates" of each individual using random numbers. Let's implement this method. We will have to iterate through individuals one by one, calculate a random number for each, and based on its value, decide whether this individual "survives" or "dies." We will start by creating a SelectionResult vector, which will reflect the results of selection. The length of this vector must be equal to the length of the Population vector. We will fill this vector with NA values, indicating missing data.
SelectionResult <- rep(NA, length(Population))
SelectionResult
## [1] NA NA NA NA NA NA NA NA NA NANow let's create a loop, and for this, we will use the function for(var in seq) expr. In this case, var is a variable that serves as the loop counter, seq is a vector from which this counter selects values, and expr is an expression that is repeated for each value of var. The expression that will be executed on each "turn" of the loop is a conditional construct if(cond) expr, where cond is the condition being checked, and expr is the command that is executed if the condition is true. For all values of the loop counter var (i.e., for all individuals in our population), we will create a random number, the probability of which is uniformly distributed between 0 and 1. Such a random number is created by the command runif(1). We will compare this random number with the corresponding value from the Mortality vector. If the random number is less than the probability of the individual dying, it should survive. In this case, we transfer the code of this individual from the Population vector to the SelectionResult vector. If the random number is less than the probability of the individual dying, this individual dies, and the NA value remains in the SelectionResult vector at the corresponding position.
for (var in 1:length(Population))
if(runif(1) < Mortality[var]) SelectionResult[var] <- Population[var]
SelectionResult
## [1] NA NA 3 NA NA NA NA 8 NA NASelection simulation is complete! In the SelectionResult vector, we see the codes of the individuals that survived, and at the positions of the individuals that died, the NA designations remain. For further model operation, NA should be removed from the results vector. This task is done by the na.omit() command.
Survivors <- na.omit(SelectionResult)
Survivors
## [1] 3 8
## attr(,"na.action")
## [1] 1 2 4 5 6 7 9 10
## attr(,"class")
## [1] "omit"You can see that when the Survivors vector, created by the na.omit() command, is displayed in the console, R reports the positions where the NA values were located. We do not need this information. It can be ignored or removed; to do this, the as.vector() command should be used to indicate that the result is needed in vector form.
Survivors <- as.vector(na.omit(SelectionResult))
Survivors
## [1] 3 8Now, in the Survivors vector, we see a list of codes of those individuals who have passed through the crucible of selection. The method described in this case is logical, understandable, and universal. Unfortunately, calculations using this method are very slow. Simulating processes in populations containing a noticeable number of individuals (say, in Type III models) will require a significant amount of time. In general, it is known that R works slowly with loops. What to do? Look for ways to achieve the same result in more economical ways that will require much less time. One possible way is to use the sample() function. 4.2.4 Features of the sample() function First of all, it is necessary to thoroughly understand the operation of the R function on which the simulation of selection will be based – the sample() function. The function sample(x, size, replace = FALSE, prob = NULL) can be used both for generating sets of random integers and for selecting a certain number of objects from a population. The first two arguments are mandatory: x is the population (vector) from which the selection is made, and size is the sample size to be obtained. The next two parameters are optional. The replace parameter determines whether an element that has already been selected can be selected again. If replace = FALSE, after a certain element is selected, the selected value is removed from those that can be selected in the future. If replace = TRUE is specified, the command will select elements from the complete population each time, regardless of whether these elements were selected before or not. The prob = NULL parameter, which should specify a vector of distribution "weights" (relative probability of selecting different elements), is interesting. This parameter can be very useful for simulating selection. In one of the variants, which we will describe later, using this parameter greatly simplifies the life of the user of simulation models. If the replace and prob parameters are not specified, R uses the values provided in the sample record of arguments: sample(x, size, replace = FALSE, prob = NULL). Let's consider how the sample function works with a simple example.
sample(1:10, 9)
## [1] 3 2 1 6 8 9 5 4 10The expression 1:10 defines a vector of integers from one to ten. Nine random objects are taken from this population. As a result, we obtained a population of 9 randomly arranged numbers out of 10 (one number is missing). Try executing the command sample(1:10, 11)! Why does R report an error? If each object is selected once, it is impossible to select 11 out of 10! But if repeated selection of numbers is allowed, the command will be executed.
sample(1:10, 11, replace=T)
## [1] 8 6 5 7 7 2 9 9 2 5 7This function can be used not only to generate random numbers but also to select objects from a given list. R has a list of 26 letters of the English alphabet, the object letters. Let's create a random sample of ten letters, which may be repeated.
sample(letters, 10, replace = T)
## [1] "b" "i" "h" "f" "s" "i" "r" "m" "k" "e"Does it look bad? Perhaps increase the frequency of vowel selection? Of course, assuming that different letters are equally probable is a very rough assumption. Data on the frequencies of English letters can be found on Wikipedia: "English alphabet." Let's create a vector let_probabs containing the letter frequencies taken from Wikipedia.
let_probabs <- c(0.0817, 0.0149, 0.0278, 0.0425, 0.1270, 0.0223, 0.0202, 0.0609, 0.0697, 0.0015, 0.0077, 0.0403, 0.0241, 0.0675, 0.0751, 0.0193, 0.0010, 0.0599, 0.0633, 0.0906, 0.0276, 0.0098, 0.0236, 0.0015, 0.0197, 0.0007)
let_probabs
## [1] 0.0817 0.0149 0.0278 0.0425 0.1270 0.0223 0.0202 0.0609 0.0697 0.0015
## [11] 0.0077 0.0403 0.0241 0.0675 0.0751 0.0193 0.0010 0.0599 0.0633 0.0906
## [21] 0.0276 0.0098 0.0236 0.0015 0.0197 0.0007
sample(letters, 10, replace = T, prob = let_probabs)
## [1] "y" "i" "e" "o" "h" "u" "t" "h" "h" "i"Now the generated sequence of letters looks better! But we can ask an important question. The let_probabs vector lists the probabilities of individual letters. The sum of probabilities of all alternatives must be 1. Is this really the case?
sum(let_probabs)
## [1] 1.0002Why is the sum of probabilities higher? We took the data on letter probabilities from Wikipedia, and they are rounded there. It turns out that rounding "up" occurred somewhat more often than rounding "down." We can draw an important conclusion. The prob parameter actually uses not probabilities (the sum of probabilities is 1), but "weights," relative probabilities! Let's try to create a file with random weights for letters – for example, such that they monotonically decrease from the beginning of the alphabet to its end.
let_weights <- 26:1
let_weights
## [1] 26 25 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2
## [26] 1
sum(let_weights)
## [1] 351
sample(letters, 10, replace = T, prob = let_weights)
## [1] "i" "l" "e" "c" "h" "a" "e" "k" "e" "d"Even though the sum of the elements in this file is 351, this method of changing the probabilities of selecting different elements is quite workable. The sample generated according to these "weights" predominantly contains letters located closer to the beginning of the alphabet. Thus, we can provide any numbers to the vector specified in the prob parameter as the source of "weights" (probability ratios). What matters is not their sum, but their ratio. Of course, the length of the source vector (from which elements are selected) and the "weights" vector must be the same. 4.2.5 Differential mortality among individuals with the same viability: using a "fate vector" In the previous case, we created a Population vector listing the codes of 10 individuals. How to simulate selection with a frequency s (let's choose the previously discussed option s = 0.3 for example) among these individuals? Let's create a Fate vector for this purpose, which will determine which individuals will die and which will survive.
s <- 0.3Let's recall the command format: sample(x, size, replace = T, prob = vect). x is the population from which elements should be selected. Since we are building a "fate vector," there can be two options: death or survival. Creating a vector containing these options is not difficult; for this, we need to combine two required values: c(NA, 1). What should be the length of the "fate vector"? The same as the vector containing the list of individuals. The number of elements in it can be determined as length(Population). In the "fate vector," values leading to death (NA) should be placed with a frequency corresponding to the probability of dying (s), and values allowing survival (1) should occur with a frequency (1-s). These two probabilities should be combined into a vector of the same length as the vector containing the selection results: c(s, 1-s). Overall, it looks like this:
Fate <- sample(c(NA, 1), length(Population), replace=T, prob=c(s, 1-s))
Fate
## [1] 1 1 1 1 1 NA 1 1 1 1The die is cast. It remains to apply it to determine the fate of each individual. And here we will use the features of R as a vectorized language. The fact is that in R, the product of vectors, if they have the same length, is a vector of the same length, where the first value is the product of the first values of the vectors, the second is the product of the second, and so on. As in the case when we simulated selection using a for loop, we get a SelectionResult vector containing NA values. As in the mentioned case, from it, we need to obtain a Survivors vector, which contains only those individuals who survived.
SelectionResult <- Population*Fate
SelectionResult
## [1] 1 2 3 4 5 NA 7 8 9 10
Survivors <- as.vector(na.omit(SelectionResult))
Survivors
## [1] 1 2 3 4 5 7 8 9 104.2.6 Differential mortality in a group consisting of different subgroups using a combined "fate vector" Let's start by creating a MixedPopul vector consisting of individuals with codes 1 and 2. During the operation of the simulation model, a similar file arises as a result of previous calculations; in this case, we will have to create it specifically. The sample() function can also be used for this.
MixedPopul <- sample(c(1, 2), 20, replace=T, prob=c(0.6, 0.4))Suppose selection with intensity s1 = 0.3 acts against individuals with code 1; and selection with intensity s2 = 0.1 acts against individuals with code 2.
s1 <- 0.3
s2 <- 0.1The MixedPopul vector has individuals with code 1 and code 2 mixed up. To determine their location, we will use the command which(x==a), which returns the ordinal elements of x that are equal to a.
MixedPopul
## [1] 1 2 1 1 2 1 2 2 2 1 2 1 1 1 2 1 1 1 2 2
which(MixedPopul==1)
## [1] 1 3 4 6 10 12 13 14 16 17 18
which(MixedPopul==2)
## [1] 2 5 7 8 9 11 15 19 20Thanks to the which() command, we have formed vectors. What is the length of these vectors?
length(which(MixedPopul==1))
## [1] 11
length(which(MixedPopul==2))
## [1] 9For individuals with code 1, we need to create a "fate vector" Fate1, where the ratio of randomly arranged NA and 1 values would depend on the value of s1. The same should be done for individuals with code 2.
Fate1 <- sample(c(NA, 1), length(which(MixedPopul==1)), replace=T, prob=c(s1, 1-s1))
Fate1
## [1] 1 NA 1 NA NA 1 1 1 NA 1 1
Fate2 <- sample(c(NA, 1), length(which(MixedPopul==2)), replace=T, prob=c(s2, 1-s2))
Fate2
## [1] 1 1 1 1 1 1 1 1 1Now our task is to combine the "fate vectors" Fate1 and Fate2 into a combined vector Fate, where at the positions where 1 is in the MixedPopul vector, values from Fate1 are located, and where 2 is, values from Fate2 are located. Let's do that.
Fate <- rep(NA, length(MixedPopul))
Fate[which(MixedPopul==1)] <- Fate1
Fate[which(MixedPopul==2)] <- Fate2
Fate
## [1] 1 1 NA 1 1 NA 1 1 1 NA 1 1 1 1 1 NA 1 1 1 1Can this sequence of commands be simplified? Of course. We can avoid creating separate Fate1 and Fate2 vectors and directly form a "fate vector" (let's call it FATE to distinguish it from the previous vector; you do remember that case matters in R!).
FATE <- rep(NA, length(MixedPopul))
FATE[which(MixedPopul==1)] <- sample(c(NA, 1), length(which(MixedPopul==1)), replace=T, prob=c(s1, 1-s1))
FATE[which(MixedPopul==2)] <- sample(c(NA, 1), length(which(MixedPopul==2)), replace=T, prob=c(s2, 1-s2))
FATE
## [1] 1 1 1 1 1 NA 1 1 1 NA 1 1 NA NA 1 1 1 NA 1 1It remains to repeat the sequence of commands that we used in previous cases!
SelectionResult <- MixedPopul*FATE
Survivors <- as.vector(na.omit(SelectionResult))
Survivors
## [1] 1 2 1 1 2 2 2 2 2 1 2 1 1 2 24.2.7 Simulating non-selective population reduction We have considered three methods for simulating differential survival. Selection reduces the population size, and reproduction increases it. In most cases, these two processes are not precisely balanced. For the model population to be stable, reproduction must add more individuals than selection removes. But how to prevent population overgrowth then? Through either non-selective reduction (reduction not affected by competition between individuals) or competitive (selective) reduction. As noted above, competitive population reduction is one manifestation of selection. We will examine it in detail in another section. In the models we will build now, we will consider the simplest case – non-selective reduction. We will perform it when transitioning from the omega-population of one cycle to the alpha-population of the next. We need to select from the omega-population of the previous cycle a number of individuals for the alpha-population of the next cycle equal to N – the population size limit. Using the mixed formula notation discussed in section "3.3 The notation system for work cycle stages and variables adopted in the course," we can describe non-selective reduction as follows: α t P = sample( ω t–1 P, N). Let's try to perform such calculations.
N <- 10
Om_Pop <- 1:12
Al_Pop <- sample(Om_Pop, N)
Al_Pop
## [1] 6 10 1 3 7 9 2 4 12 11The provided script will not work correctly in all cases. The situation where N > length(Om_Pop) should be anticipated – could selection have thinned the population significantly? In that case, this method will not work: try executing the content of the previous chunk (programmatic insertion in the text), where N <- 15 is specified. What to do then? One possible solution is to use the conditional construct if... else.... Let's examine it in more detail: if(cond) cons.expr else alt.expr, where cond is the logical condition being checked, cons.expr is the action performed if the condition is true, and alt.expr is the action performed when the condition is false.
if (N >= length(Om_Pop)) Al_Pop <- Om_Pop else Al_Pop <- sample(Om_Pop, N)Let's check how this works. We just created the Om_Pop vector with 12 elements.
N <- 15
if (N >= length(Om_Pop)) Al_Pop <- Om_Pop else Al_Pop <- sample(Om_Pop, N)
Al_Pop
## [1] 1 2 3 4 5 6 7 8 9 10 11 12
N <- 12
if (N >= length(Om_Pop)) Al_Pop <- Om_Pop else Al_Pop <- sample(Om_Pop, N)
Al_Pop
## [1] 1 2 3 4 5 6 7 8 9 10 11 12
N <- 10
if (N >= length(Om_Pop)) Al_Pop <- Om_Pop else Al_Pop <- sample(Om_Pop, N)
Al_Pop
## [1] 12 4 9 5 6 8 7 11 10 14.3 Implementation of the conceptual model of monogenic selection in R 4.3.1 Calculation algorithm Now that we have established how to simulate differential mortality and non-selective reduction, we can implement what is described in section 4.1.3, "Sequence of calculated quantities in the model's work cycle." When compiling the algorithm, we will use the mixed formula notation. ENTRANCE: Initial state of the system – initial state of the modeled system initial population sizes of both forms: α 1 N A and α 1 N a; Parameters – parameters population size: N; Malthus parameter: r (we will choose exclusively natural values); selection coefficient s, where s ≈ 1 – ( t+1 N A / t N A ) / ( t+1 N a / t N a ). Experimental conditions – model experiment conditions number of model cycles: cycles; (if necessary, for Type II model): number of iterations (runs of the model with the same initial conditions) iterat; Changeable parameters combinations – combinations of changeable parameters (if necessary, for Type III model): sets of initial parameter values tested in possible combinations. CALCULATIONS: Alpha population (breeders): If t = 1, the population sizes are given by the initial parameters, **** α 1 P = α 1 N A + α 1 N a **; If t > 1, the alpha population is formed during non-selective reduction of the omega population: if (N >= length( ω t-1 P)) α^ t P <- ω t-1 P else α^ t P <- sample( ω t-1 P, N). Beta population (offspring): β t P = r × α t P. Omega population (offspring thinned by selection): t Fate <- rep(1, length( β t P)) t Fate[which( β t P==A)] <- sample(c(NA, 1), length(which( β t P==A)), replace=T, prob=c(s, 1-s)) ω t P <- as.vector(na.omit( β t P * t Fate)) In fact, the model is already created; everything else is the implementation of the conceived plan. 4.3.2 Creating a Type I model: initial parameters and experiment conditions Open RStudio, create a new file (R script), and start writing the model! First, set the informational input to the model: initial state, parameters, and experiment conditions. Where to get the values to specify in the model? From anywhere, even "out of thin air"! To write the script and observe its operation, some values need to be set. How meaningful they are, how they relate to the goal of the model experiment, will be decided later. Initially, it is better to choose small values to make it easier to see intermediate calculation results in the console.
# ДОБІР ЗА ОДНІЄЮ ОЗНАКОЮ В ПОПУЛЯЦІЇ МОНОЦИКЛІЧНИХ КЛОНАЛЬНИХ ОРГАНІЗМІВ З СИНХРОННИМ РОЗМНОЖЕННЯМ (Модель І типу)
# Популяція складається з двох форм: A (код 1) та a (код 2). Проти A діє добір s
# Робочий цикл: альфа-популяція → розмноження (r) → бета → скорочення (N) → омега → альфа
# Initial script commands — початкові команди скрипта
# setwd("~/data/SexOnR") # Робоча директорія (лише для Д.Ш.!!!)
rm(list = ls()) # Очищення раніше збережених об'єктів в Environment
set.seed(123456) # Забезпечення повторюваності у виборі псевдовипадкових чисел
# I. ENTRANCE — ВХІД:
# Initial state of the system — початковий стан модельованої системи
al_1_N_a <- 1 # Початкова чисельність a-особин
# Parameters — параметри
N <- 20 # Альфа-чисельність модельної популяції
al_1_N_A <- N - al_1_N_a # Початкова чисельність A-особин (усі крім a-особин)
N <- 10 # Альфа-чисельність модельної популяції
r <- 2 # Параметр Мальтуса (кількість потомків на одного плідника на один цикл роботи моделі)
s <- 0.1 # Коефіцієнт добору на користь a-особин
# Experimental conditions — умови експерименту з моделлю
cycles <- 6 # Кількість робочих циклів моделі Initial parameters can be set as specific numbers or as certain dependencies (formulas), as we did for α 1 N A. If we are simulating selection in favor of form a, it is logical to start with a case where the proportion of a in the model population (its "concentration") is minimal. This line (which uses the N parameter) should be placed after the N parameter has been defined; this is why it "jumps" from the "Initial state of the system" section of the script to "Parameters." When creating a model, it is advisable to set relatively small values for parameters and cycle counts. When writing the model step-by-step (as will be shown later), it is useful to check what comes out at each step. In this case, it is convenient for the objects that will be output to the console not to be too large. 4.3.3 Creating objects that the model will work with Planning the objects to be created in the model should start "from the end," from the form in which the simulation result will be obtained. In what form will the results of the model's work be output? How should the data be organized? Most likely, the main result of a certain iteration is the dynamics of the "concentration" of one of the forms over time, i.e., changes in t C a depending on t. For a single iteration, this dynamic can be stored in one vector; for a series of consecutive iterations, a matrix will be needed. Let's start with the variant that we will complicate further (we chose it when setting the number of iterations greater than one in the initial conditions).
# II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
Result <- rep(NA, cycles+1) # Створення вектора для запису підсумків
Result
## [1] NA NA NA NA NA NA NABefore starting the model's work cycle, it is advisable to create the objects that will be used in it. R, in general, is a language not very adapted to working with loops. This language processes loops where new objects are created most slowly. Therefore, before starting the model's work cycle, it is advisable to create the objects it will work with. Overwriting values in a ready-made object is a simpler task than creating a new object. However, for some objects, an exception will have to be made, and they will have to be created inside the loop. The fact is that the length of these objects will be different for different cycles. But at this stage, creating the model, we can make it suboptimal. We will create all objects, except for the vector for the final result, inside the loop. Later, when the model is running, it can be optimized (especially when converting it to a Type III model). However, it makes sense to create the vector for the initial composition at each cycle (Al_Pop) immediately. Its length, according to the conditions provided by the algorithm, should be equal to N. First, let's make this file empty, i.e., filled with NA values.
Al_Pop <- rep(NA, N) # Створення файлу для початкового складуThe alpha composition of the population in the first cycle is set differently than in all subsequent ones. Therefore, its initial definition can be moved forward, before the start of the work cycle. Since this quantity is calculated separately, its recording in the results file should also be done separately.
# III. CALCULATIONS — РОЗРАХУНКИ:
# Initial composition creation — утворення початкового складу
Al_Pop[1:al_1_N_a] <- 2 # Перенесення до початкового складу a-особин
Al_Pop[(al_1_N_a+1):(al_1_N_a+al_1_N_A)] <- 1 # Перенесення до початкового складу A-особин
Al_Pop
## [1] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1As you can see, we denoted form A as 1 and form a as 2. The method used only works if 1 N a ≥ 1. Should you always output the created objects to the console? Of course not. But right now, as we create the script step by step, it is very useful: you can immediately see if you are getting what you need. When writing a script, it is very easy to make mistakes; constant feedback helps to find and correct them quickly. If the alpha population of the first cycle has been created, its characteristic can be recorded in the results file.
Result[1] <- length(which(Al_Pop==2)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==2))) # Запис до вектора з підсумками першої позиції: початкового складу на початку роботи моделі
Result
## [1] 0.05 NA NA NA NA NA NA4.3.4 Stages of the model's working cycle Objects are created, the working cycle can begin. As we create it, we want to get feedback: to check if we are getting what we planned. How to do this when the cycle is not yet fully written? Plan a place where the cycle will start, but temporarily close it with a '#'. The cycle should iterate through the values of the variable t. In this case, let's mark it manually! When the cycle is fully written and works, we will remove these temporary constructs.
# Main work cycle — основний робочий цикл
# for (t in 1:cycles) { # Робочий цикл моделі (його початок тимчасово вимкнено)
t <- 1 # після запуску циклу цю строку можна буде заблокувати або прибрати
# (тут будемо покроково створювати цикл)
# } # Кінець робочого циклу (тимчасово вимкнено)When it comes to the first cycle t=1, the alpha population has already been created before the start of the model's working cycle. What about all other cases? Paradoxically, it can be created at the end of the working cycle. In this case, in any case, the working cycle will begin with the calculation of the beta population based on the already existing alpha population. Let's start. We will create a vector for offspring and see what happens.
Be_Pop <- rep(Al_Pop, r) # Заповнення вектора для потомства
Be_Pop
## [1] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [39] 1 1The next stage is selection. We need to create a 'fate vector' Fate. Let's fill it with ones (i.e., values that will ensure the survival of the corresponding individuals from the beta population). NA, which will denote death, can be inserted later specifically for those individuals who are destined to die. After that, we will use the calculation method discussed in section 4.2.6 'Differential mortality in a group consisting of different subgroups using a combined 'fate vector''. For those positions in the 'fate vector' where A-individuals are located, we will place NA and 1 symbols in the ratio determined by the selection coefficient s.
Fate <- rep(1, length(Be_Pop)) # Створення «вектору долі»
Fate[which(Be_Pop==1)] <- sample(c(NA, 1), length(which(Be_Pop==1)), replace=T, prob=c(s, 1-s)) # Деякі з A-особин можуть загинути
Fate
## [1] 1 1 1 1 1 1 1 1 1 NA 1 1 1 NA 1 NA 1 1 1 1 1 1 1 1 1
## [26] 1 1 1 1 1 1 1 1 1 1 1 1 1 NA 1Now we can calculate the omega population.
Om_Pop <- as.vector(na.omit(Be_Pop * Fate)) # Визначення, які особини виживають, а які гинуть
Om_Pop
## [1] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1All that remains is to carry out a non-selective reduction in the size of the model population to N.
if (N >= length(Om_Pop)) Al_Pop <- Om_Pop else Al_Pop <- sample(Om_Pop, N) # Створення альфа-популяції для наступного цикла
Al_Pop
## [1] 1 1 1 1 1 1 1 1 1 1Are all components of the model's working cycle defined? No. In the previous step of cycle t, we obtained the alpha composition of the model population for cycle t+1. It should be written into the Result matrix. But before that, two more rows should be added, one of which will work if only one form remains in the model population. If there were a-individuals in the alpha composition of the previous cycle, they should have entered the omega composition. But if there are few of them, due to simple probability, they may not enter the alpha count of the next cycle. Is it worth carrying out further calculations in this case? The operation of type I model can be stopped, and in type II model, you can proceed to the next iteration. If only individuals of one form remain in the model population, the composition of this population will no longer change. The work (in any case, in this iteration) should be stopped and it should be recorded that it ended with the disappearance of a-individuals. Where to record? In the last row of the Result matrix! Similarly, the possibility of the disappearance of A-individuals from the model population should be foreseen; they could have disappeared during the transition from the beta population to the omega population or further, during the transition to the alpha population. To foresee these possibilities, we will use if (condition) { ...; ... } conditions. The condition to be checked is placed in parentheses. If this condition is true, the subsequent commands are executed. A series of commands to be executed sequentially can be grouped in curly braces and either written on consecutive lines or separated by semicolons. The break function interrupts the execution of the for loop, i.e., ensures the interruption of iterating through consecutive working cycles. When running the next fragment of R code, R Markdown does not 'understand' which loop the break commands belong to, as the model's working cycle is not currently running (we haven't turned it on yet). Let's close the lines with these commands with '#'.
# «Запобіжники» на випадок зникнення однієї з форм:
# if (length(which(Al_Pop==1))==0) {Result[cycles+1] <- 1; break} # У разі перемоги a-особин
# if (length(which(Al_Pop==2))==0) {Result[cycles+1] <- 0; break} # У разі перемоги A-особин
Result[t+1] <- length(which(Al_Pop==2)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==2))) # Запис наступного рядка до вектора з результатами роботи
# } # Кінець робочого циклуAlthough the model's working cycle has not yet been launched, corresponding entries should have appeared in the Result vector at the first and second positions. 4.3.5 Launching the model's working cycle Do you remember where we marked the start of the working cycle in our script? Now it needs to be launched. Let's present the entire model, which does not yet include visualization of the results. We will use the following input conditions, which will make it more interesting to observe the model.
# ДОБІР ЗА ОДНІЄЮ ОЗНАКОЮ В ПОПУЛЯЦІЇ МОНОЦИКЛІЧНИХ КЛОНАЛЬНИХ ОРГАНІЗМІВ З СИНХРОННИМ РОЗМНОЖЕННЯМ (Модель І типу)
# Популяція складається з двох форм: A (код 1) та a (код 2). Проти A діє добір s
# Робочий цикл: альфа-популяція → розмноження (r) → бета → скорочення (N) → омега → альфа
# Initial script commands — початкові команди скрипта
# setwd("~/data/SexOnR") # Робоча директорія (лише для Д.Ш.!!!)
rm(list = ls()) # Очищення раніше збережених об'єктів в Environment
set.seed(123456) # Забезпечення повторюваності у виборі псевдовипадкових чисел
# I. ENTRANCE — ВХІД:
# Initial state of the system — початковий стан модельованої системи
al_1_N_a <- 1 # Початкова чисельність a-особин
# Parameters — параметри
N <- 25 # Альфа-чисельність модельної популяції
al_1_N_A <- N - al_1_N_a # Початкова чисельність A-особин (усі крім a-особин)
r <- 2 # Параметр Мальтуса (кількість потомків на одного плідника на один цикл роботи моделі)
s <- 0.2 # Коефіцієнт добору на користь a-особин
# Experimental conditions — умови експерименту з моделлю
cycles <- 20 # Кількість робочих циклів моделі
# II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
Result <- rep(NA, cycles+1) # Створення вектора для запису підсумків
Al_Pop <- rep(NA, N) # Створення файлу для початкового складу
# III. CALCULATIONS — РОЗРАХУНКИ:
# Initial composition creation — утворення початкового складу
Al_Pop[1:al_1_N_a] <- 2 # Перенесення до початкового складу a-особин
Al_Pop[(al_1_N_a+1):(al_1_N_a+al_1_N_A)] <- 1 # Перенесення до початкового складу A-особин
Result[1] <- length(which(Al_Pop==2)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==2))) # Запис до вектора з підсумками першої позиції: початкового складу на початку роботи моделі
# Main work cycle — основний робочий цикл
for (t in 1:cycles) { # Робочий цикл моделі
Be_Pop <- rep(Al_Pop, r) # Заповнення вектора для потомства
Fate <- rep(1, length(Be_Pop)) # Створення «вектору долі»
Fate[which(Be_Pop==1)] <- sample(c(NA, 1), length(which(Be_Pop==1)), replace=T, prob=c(s, 1-s)) # Деякі з A-особин можуть загинути
Om_Pop <- as.vector(na.omit(Be_Pop * Fate)) # Визначення, які особини виживають, а які гинуть
if (N >= length(Om_Pop)) Al_Pop <- Om_Pop else Al_Pop <- sample(Om_Pop, N) # Створення альфа-популяції для наступного цикла
# «Запобіжники» на випадок зникнення однієї з форм:
if (length(which(Al_Pop==1))==0) {Result[cycles+1] <- 1; break} # У разі перемоги a-особин
if (length(which(Al_Pop==2))==0) {Result[cycles+1] <- 0; break} # У разі перемоги A-особин
Result[t+1] <- length(which(Al_Pop==2)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==2))) # Запис наступного рядка до вектора з результатами роботи
} # Кінець робочого циклу
# IV. FINISHING — ЗАВЕРШЕННЯ:
# Results viewing — перегляд підсумків
Result # Перегляд вектора з динамікою концентрації a-особин
## [1] 0.04 0.04 0.04 0.04 0.08 0.08 0.04 0.04 0.08 0.12 0.16 0.28 0.40 0.48 0.56
## [16] 0.60 0.72 0.84 0.88 0.92 1.00It works!!! All that's left is to add the graph.
# Results visualization — візуалізація підсумків
plot(Result, type="l", lty=1, col="red") # ГрафікFig. 4.3.5.1 Visualization of the dynamics of the proportion of a-individuals in the Type I model The model (in its minimal configuration) has started working. But this is far from everything. 4.3.6 Example of Type I model operation Let's run the program in the variant shown below (increase the population size and the number of cycles to see a larger picture than during model debugging). Let's add an 'Initial script commands' block to the script. Note: the set.seed() command is active, meaning that each time the program is executed, it will lead to the same result. The last line is a command to plot the results of the first iteration on the graph.
# ДОБІР ЗА ОДНІЄЮ ОЗНАКОЮ В ПОПУЛЯЦІЇ МОНОЦИКЛІЧНИХ КЛОНАЛЬНИХ ОРГАНІЗМІВ З СИНХРОННИМ РОЗМНОЖЕННЯМ (Модель І типу)
# Популяція складається з двох форм: A (код 1) та a (код 2). Проти A діє добір s
# Робочий цикл: альфа-популяція → розмноження (r) → бета → скорочення (N) → омега → альфа
# Initial script commands — початкові команди скрипта
# setwd("~/data/SexOnR") # Робоча директорія (лише для Д.Ш.!!!)
rm(list = ls()) # Очищення раніше збережених об'єктів в Environment
set.seed(123456) # Забезпечення повторюваності у виборі псевдовипадкових чисел
# I. ENTRANCE — ВХІД:
# Initial state of the system — початковий стан модельованої системи
al_1_N_a <- 1 # Початкова чисельність a-особин
# Parameters — параметри
N <- 100 # Альфа-чисельність модельної популяції
al_1_N_A <- N - al_1_N_a # Початкова чисельність A-особин (усі крім a-особин)
r <- 2 # Параметр Мальтуса (кількість потомків на одного плідника на один цикл роботи моделі)
s <- 0.02 # Коефіцієнт добору на користь a-особин
# Experimental conditions — умови експерименту з моделлю
cycles <- 100 # Кількість робочих циклів моделі
# II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
Result <- rep(NA, cycles+1) # Створення вектора для запису підсумків
Al_Pop <- rep(NA, N) # Створення файлу для початкового складу
# III. CALCULATIONS — РОЗРАХУНКИ:
# Initial composition creation — утворення початкового складу
Al_Pop[1:al_1_N_a] <- 2 # Перенесення до початкового складу a-особин
Al_Pop[(al_1_N_a+1):(al_1_N_a+al_1_N_A)] <- 1 # Перенесення до початкового складу A-особин
Result[1] <- length(which(Al_Pop==2)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==2))) # Запис до вектора з підсумками першої позиції: початкового складу на початку роботи моделі
# Main work cycle — основний робочий цикл
for (t in 1:cycles) { # Робочий цикл моделі
Be_Pop <- rep(Al_Pop, r) # Заповнення вектора для потомства
Fate <- rep(1, length(Be_Pop)) # Створення «вектору долі»
Fate[which(Be_Pop==1)] <- sample(c(NA, 1), length(which(Be_Pop==1)), replace=T, prob=c(s, 1-s)) # Деякі з A-особин можуть загинути
Om_Pop <- as.vector(na.omit(Be_Pop * Fate)) # Визначення, які особини виживають, а які гинуть
if (N >= length(Om_Pop)) Al_Pop <- Om_Pop else Al_Pop <- sample(Om_Pop, N) # Створення альфа-популяції для наступного цикла
# «Запобіжники» на випадок зникнення однієї з форм:
if (length(which(Al_Pop==1))==0) {Result[cycles+1] <- 1; break} # У разі перемоги a-особин
if (length(which(Al_Pop==2))==0) {Result[cycles+1] <- 0; break} # У разі перемоги A-особин
Result[t+1] <- length(which(Al_Pop==2)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==2))) # Запис наступного рядка до вектора з результатами роботи
} # Кінець робочого циклу
# IV. FINISHING — ЗАВЕРШЕННЯ:
# Results viewing — перегляд підсумків
Result # Перегляд вектора з динамікою концентрації a-особин
## [1] 0.01 0.02 0.03 0.03 0.04 0.05 0.04 0.05 0.08 0.11 0.16 0.17 0.22 0.26 0.25
## [16] 0.23 0.21 0.21 0.22 0.24 0.22 0.20 0.20 0.18 0.14 0.12 0.16 0.18 0.20 0.17
## [31] 0.16 0.14 0.16 0.18 0.15 0.14 0.15 0.12 0.14 0.17 0.15 0.16 0.15 0.14 0.11
## [46] 0.08 0.08 0.08 0.07 0.02 0.04 0.05 0.06 0.06 0.05 0.07 0.08 0.11 0.09 0.07
## [61] 0.05 0.02 NA NA NA NA NA NA NA NA NA NA NA NA NA
## [76] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [91] NA NA NA NA NA NA NA NA NA NA 0.00
# Results visualization — візуалізація підсумків
plot(Result, type="l", lty=1, col="red", ylim=c(0, 1),
main = "The part of the form a in the population a and A",
xlab = "Simulation cycles", ylab = "The part of the form a") # ГрафікFig. 4.3.6.1 Another example of the dynamics of the proportion of a-individuals in the Type I model. Graph and axis titles added What do we have? The proportion ('concentration') of a-individuals, shown by the red line, increased rather than decreased for a period of time, and then decreased rather than increased. It ended with all a-individuals disappearing from the model population. Such dynamics of the studied variable can be called a random walk: it is determined mainly by random circumstances, not by the action of selection. So, under the conditions used in the model, can selection change the composition of the model population? It can, but not always... 4.4 Monogenic selection: Type II and III models 4.4.1 Rebuilding the Type I model into a Type II model As we noted in Chapter 3, Type I models perform one iteration (one pass of the model under certain conditions), Type II models perform a group of iterations under the same conditions, and Type III models perform groups of iterations corresponding to different combinations of studied parameters. In subsection 4.3, we created a Type I model. What is needed to transform it into a Type II model? It is necessary to add tools that ensure the performance of a group of iterations under the same conditions. This is not difficult. We will 'stretch' the iteration loop over the model's working cycle. What changes in the Type I model are needed for this? First of all, if we plan to get an object that will contain all the results (changes in the 'concentration' of individuals favored by selection), it will be inconvenient to deal with a vector. It is more convenient to make a matrix where each column corresponds to a specific iteration. If the Result object becomes a matrix, references to it must be two-dimensional. It is not enough for us to indicate which cycle (t) the value we are dealing with belongs to; we also need to indicate which iteration this value belongs to: Result[t, i]. In all other respects, our Type II model does not differ from the Type I model.
# ДОБІР ЗА ОДНІЄЮ ОЗНАКОЮ В ПОПУЛЯЦІЇ МОНОЦИКЛІЧНИХ КЛОНАЛЬНИХ ОРГАНІЗМІВ З СИНХРОННИМ РОЗМНОЖЕННЯМ (Модель ІІ типу)
# Групи ітерацій. Популяція складається з двох форм: A (код 1) та a (код 2). Проти A діє добір s
# Робочий цикл: альфа-популяція → розмноження (r) → бета → скорочення (N) → омега → альфа
# Initial script commands — початкові команди скрипта
# setwd("~/data/SexOnR") # Робоча директорія (лише для Д.Ш.!!!)
rm(list = ls()) # Очищення раніше збережених об'єктів в Environment
set.seed(123456) # Забезпечення повторюваності у виборі псевдовипадкових чисел
# I. ENTRANCE — ВХІД:
# Initial state of the system — початковий стан модельованої системи
al_1_N_a <- 3 # Початкова чисельність a-особин
# Parameters — параметри
N <- 100 # Альфа-чисельність модельної популяції
al_1_N_A <- N - al_1_N_a # Початкова чисельність A-особин (усі крім a-особин)
r <- 2 # Параметр Мальтуса (кількість потомків на одного плідника на один цикл роботи моделі)
s <- 0.1 # Коефіцієнт добору на користь a-особин
# Experimental conditions — умови експерименту з моделлю
cycles <- 50 # Кількість робочих циклів моделі
iterat <- 3 # Кількість ітерацій
# II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
Result <- matrix(NA, nrow = cycles+1, ncol = iterat) # У разі використання ітерацій для запису результатів потрібна матриця
Al_Pop <- rep(NA, N) # Створення файлу для початкового складу
# III. CALCULATIONS — РОЗРАХУНКИ:
# Higher level cycles running — запуск циклів вищого рівня
for (i in 1:iterat) { # Початок ітерації
# Initial composition creation — утворення початкового складу
Al_Pop[1:al_1_N_a] <- 2 # Перенесення до початкового складу a-особин
Al_Pop[(al_1_N_a+1):(al_1_N_a+al_1_N_A)] <- 1 # Перенесення до початкового складу A-особин
Result[1, i] <- length(which(Al_Pop==2)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==2))) # Запис до вектора з підсумками першої позиції: початкового складу на початку роботи моделі
# Main work cycle — основний робочий цикл
for (t in 1:cycles) { # Робочий цикл моделі
Be_Pop <- rep(Al_Pop, r) # Заповнення вектора для потомства
Fate <- rep(1, length(Be_Pop)) # Створення «вектору долі»
Fate[which(Be_Pop==1)] <- sample(c(NA, 1), length(which(Be_Pop==1)), replace=T, prob=c(s, 1-s)) # Деякі з A-особин можуть загинути
Om_Pop <- as.vector(na.omit(Be_Pop * Fate)) # Визначення, які особини виживають, а які гинуть
if (N >= length(Om_Pop)) Al_Pop <- Om_Pop else Al_Pop <- sample(Om_Pop, N) # Створення альфа-популяції для наступного циклу
# «Запобіжники» на випадок зникнення однієї з форм:
if (length(which(Al_Pop==1))==0) {Result[cycles+1, i] <- 2; break} # У разі перемоги a-особин
if (length(which(Al_Pop==2))==0) {Result[cycles+1, i] <- 1; break} # У разі перемоги A-особин
Result[t+1, i] <- length(which(Al_Pop==2)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==2))) # Запис результата у матрицю
} # Кінець робочого циклу
# Ending of higher level cycles — завершення циклів вищого рівня
} # Кінець ітерації
# IV. FINISHING — ЗАВЕРШЕННЯ:
# Results viewing — перегляд підсумків
tail(Result) # Перегляд матриці з підсумками
## [,1] [,2] [,3]
## [46,] NA 0.16 NA
## [47,] NA 0.22 NA
## [48,] NA 0.23 NA
## [49,] NA 0.26 NA
## [50,] NA 0.25 NA
## [51,] 2 0.31 1
# Results visualization — візуалізація підсумків
plot(Result[, 1], type="l", lty=1, col="red", ylim=c(0, 1),
main = "The part of the form a in the population a and A\nin 3 different iterations",
xlab="Simulation cycles", ylab="The part of the form a")
lines(Result[, 2], type="l", lty=1, col="blue")
lines(Result[, 3], type="l", lty=1, col="green")Fig. 4.4.1.1 Results of three iterations in the Type II model We used the command set.seed(123456) and are now confident that each time the R Markdown file is executed, the result will be the same. As can be seen, out of three iterations, in the first, a-individuals favored by selection won, in the third, they disappeared, and in the second, the competition between a- and A-individuals continues. Experiment with different values of the selection coefficient, population size, and initial numbers of a-individuals! The result of the Type II model is no longer the result of a particular iteration (such results can differ significantly due to random reasons), but the probability distribution of results under certain experimental conditions. How to compare probability distributions of results? One way is as we just did. Change the parameters whose influence we want to establish and compare the results of different groups of iterations. But we are working with a programming language that has extensive capabilities! Let's create a model for comparing probability distributions using R's tools. 4.4.2 Type III model We are moving on to using the simulation model in a way that can quickly exhaust the computational resources of the computer it runs on. However, it is precisely this use of simulation models that makes them truly interesting. We have confirmed that the strength of selection and, it seems, the population size significantly affect the simulation result (in small populations, the influence of random events on their fate increases). How to investigate the influence of these parameters on the probability distribution of results? Add two more levels of loops to the model! 'On top' of the iteration loop (within which the model's working cycle operates), you can add one or two more (more - it will be difficult to analyze) loops that iterate through the values of the studied parameters. This can be done as follows.
# ДОБІР ЗА ОДНІЄЮ ОЗНАКОЮ В ПОПУЛЯЦІЇ МОНОЦИКЛІЧНИХ КЛОНАЛЬНИХ ОРГАНІЗМІВ З СИНХРОННИМ РОЗМНОЖЕННЯМ (Модель ІІІ типу)
# I. ENTRANCE — ВХІД:
# Initial script commands — початкові команди скрипта
# install.packages("plot3D")
library(plot3D)
# setwd("~/data/SexOnR") # Робоча директорія (лише для Д.Ш.!!!)
rm(list = ls()) # Очищення раніше збережених об'єктів в Environment
set.seed(123456) # Забезпечення повторюваності у виборі псевдовипадкових чисел
# Initial state of the system — початковий стан модельованої системи
al_1_N_a <- 3 # Початкова чисельність a-особин
# Parameters — параметри
# N <- 100 # Альфа-чисельність модельної популяції; ЗМІНЮВАНИЙ ПАРАМЕТР!
r <- 2 # Параметр Мальтуса (кількість потомків на одного плідника на один цикл роботи моделі)
# s <- 0.01 # Коефіцієнт добору на користь a-особин; ЗМІНЮВАНИЙ ПАРАМЕТР!
# Experimental conditions — умови експерименту з моделлю
cycles <- 200 # Кількість робочих циклів моделі
iterat <- 20 # Кількість ітерацій
# Changeable parameters combinations — комбінації змінюваних параметрів
# Створення векторів зі значеннями змінних параметрів, які будуть перебиратися у ході роботи:
vector_N <- c(100, 200, 400, 700, 1000, 1800, 3600) # Варіанти численності
vector_s <- c(0.01, 0.02, 0.04, 0.08, 0.15, 0.3) # Варіанти добору
# II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Parameters combinations mechanism — механізм комбінації змінюваних параметрів
Outcome <- matrix(NA, nrow = length(vector_N), ncol = length(vector_s), dimnames = list(vector_N, vector_s))
# Цикли, що перебирають досліджувані сполучення параметрів:
for (a in 1:length(vector_N)){ # Цикл, що змінює чисельність популяції
N <- vector_N[a] # Змінюване значення численності популяції
for (b in 1:length(vector_s)){ # Цикл, що змінює силу (коефіцієнт) добору
s <- vector_s[b] # Змінюване значення сили добору
al_1_N_A <- N - al_1_N_a # Початкова чисельність A-особин: ПІСЛЯ ВИЗНАЧЕННЯ N!!!
# Objects creation — створення об’єктів
Result <- matrix(NA, nrow = cycles+1, ncol = iterat) # У разі використання ітерацій для запису результатів потрібна матриця
Al_Pop <- rep(NA, N) # Створення файлу для початкового складу
# III. CALCULATIONS — РОЗРАХУНКИ:
# Higher level cycles running — запуск циклів вищого рівня
for (i in 1:iterat) { # Початок ітерації
# Initial composition creation — утворення початкового складу
Al_Pop[1:al_1_N_a] <- 2 # Перенесення до початкового складу a-особин
Al_Pop[(al_1_N_a+1):(al_1_N_a+al_1_N_A)] <- 1 # Перенесення до початкового складу A-особин
Result[1, i] <- length(which(Al_Pop==2)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==2))) # Запис до вектора з підсумками першої позиції: початкового складу на початку роботи моделі
# Main work cycle — основний робочий цикл
for (t in 1:cycles) { # Робочий цикл моделі
Be_Pop <- rep(Al_Pop, r) # Заповнення вектора для потомства
Fate <- rep(1, length(Be_Pop)) # Створення «вектору долі»
Fate[which(Be_Pop==1)] <- sample(c(NA, 1), length(which(Be_Pop==1)), replace=T, prob=c(s, 1-s)) # Деякі з A-особин можуть загинути
Om_Pop <- as.vector(na.omit(Be_Pop * Fate)) # Визначення, які особини виживають, а які гинуть
if (N >= length(Om_Pop)) Al_Pop <- Om_Pop else Al_Pop <- sample(Om_Pop, N) # Створення альфа-популяції для наступного циклу
# «Запобіжники» на випадок зникнення однієї з форм:
if (length(which(Al_Pop==1))==0) {Result[cycles+1, i] <- 1; break} # У разі перемоги a-особин
if (length(which(Al_Pop==2))==0) {Result[cycles+1, i] <- 0; break} # У разі перемоги A-особин
Result[t+1, i] <- length(which(Al_Pop==2)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==2))) # Запис результату у матрицю
} # Кінець робочого циклу
# Ending of higher level cycles — завершення циклів вищого рівня
} # Кінець ітерації
Outcome[a, b] <- mean(Result[cycles+1, ]) # Запис певного значення у матрицю з комбінаціями параметрів
} # Кінець циклу, що змінює силу (коефіцієнт) добору
} # Кінець циклу, що змінює чисельність популяції
# IV. FINISHING — ЗАВЕРШЕННЯ:
# Creating objects for results integration — створення об’єктів, що інтегрують результати
# save(Outcome, file = "Outcome.RData")
# Results viewing — перегляд підсумків
Outcome # Перегляд частки перемог a-особин при різних комбінаціях параметрів
## 0.01 0.02 0.04 0.08 0.15 0.3
## 100 0.113500000 0.15700000 0.3305000 0.7500000 0.90 1
## 200 0.107750000 0.13875000 0.3930000 0.8000000 0.80 1
## 400 0.029375000 0.25025000 0.4151250 0.6500000 0.90 1
## 700 0.028857143 0.18464286 0.3869286 0.6500000 0.85 1
## 1000 0.013400000 0.04855000 0.3655500 0.6500000 1.00 1
## 1800 0.008222222 0.07180556 0.3478611 0.6500000 0.90 1
## 3600 0.006041667 0.03318056 0.2882083 0.5498889 1.00 1
# Results visualization — візуалізація підсумківWhat has changed (besides the name, of course)? A few initial lines have been added. The plot3D library is activated: it will be needed. Two vectors have been added, listing the population size and selection values whose combinations will be studied: vector_N and vector_s. A matrix Outcome has been created, where one cell is provided for each combination of the two changing parameters. Two nested loops are placed on top of the iteration loops (into which the model's working cycles are nested). The first iterates through N values, the second through s values. For each combination of these values, the model performs the specified number of iterations. Note another change: the line defining one of the initial parameters, al_1_N_A <- N - al_1_N_a, had to be moved inside the loop that changes the value of N. If this line is left where it was, it can cause an error in the model's operation. After calculating a group of iterations with a specific parameter combination, its result is written to the Outcome matrix. After the calculations are finished, it is interesting to review this matrix (and what came out can be understood simply from the numbers in the matrix). The Outcome matrix becomes even clearer thanks to visualization using the plot3D package (which was activated at the beginning of the script). The command used for visualization can be ignored for now.
# Results visualization — візуалізація підсумків
persp3D(z = Outcome, x = vector_N, y = vector_s, theta = 330, phi = 30, # Побудова графіку
zlab = "\na-victory part", xlab = "N", ylab = "\ns",
expand = 0.6, facets = FALSE, scale = TRUE, ticktype = "detailed",nticks = 7, cex.axis=0.9,
clab = "a-victory", colkey = list(side = 2, length = 0.5))Fig. 4.3.6.1 Results of the Type III model operation: probability of a-individual wins for different selection intensities and population sizes We see a three-dimensional diagram. On the 'horizontal' plane are the parameters whose combinations we studied. On the 'vertical' axis is the proportion of wins for a-individuals. The color scale on the left helps to understand which value of the studied variable is characteristic for each parameter combination. We will not go into the details of the command for creating the graph now, but we will emphasize that the expression \n in the titles and labels ensures that the text is moved to the next line. Thanks to the operation of the used model, it became clear that the probability of the favored form winning is influenced almost exclusively by the strength of selection. In small populations, the role of randomness increases, but the ratio of 'wins' to 'losses' changes relatively little (compared to large populations).