11.1.1 Задача для імітаційного моделювання

Побудову першої моделі, що стосується теми курсу, ми будемо описувати та пояснювати дуже детально. Ми назвали нашу модель “добір за однією спадковою ознакою у популяції моноциклічних клональних організмів з синхронним розмноженням”. Головне слово тут — добір, а усе інше — уточнення, що стосуються особливості популяції, в якій цей добір спостерігається. Прокоментуємо ці поняття:

• добір — диференціальне (тобто таке, ймовірність якого відрізняється для представників порівнюваних груп) відтворення; добір може бути наслідком як диференціального виживання, так і/або диференціального розмноження;

• популяція — у даному випадку — сукупність особин, що спільно відтворюються та конкурують за спільний життєвий простір (обмеження якого моделюються шляхом встановлення певної чисельності популяції);

• …за однією спадковою ознакою — мова йде про добір між організмами, що належать до однієї з двох форм; важливою особливістю популяції, де може спостерігатися такий добір, є те, що ніякі інші ознаки досліджуваних організмів не впливають на їхню ймовірність вижити або розмножитися (звісно, у природних популяціях такого не буває: це абстракція); причина такої відмінності може бути, наприклад, у виникненні у чистій лінії (сукупності генетично ідентичних особин, скоріше за усе — потомків однієї особини) мутації, що стосується одного гена, та конкуренції між двома алелями: початковим та зміненим;

• моногенні ознаки — ознаки організмів, стан яких визначається одним геном; насправді, у високоорганізованих організмів таких ознак небагато, адже більшість генів водночас впливає на велику кількість ознак;

• гаплоїдні — організми, клітини тіла котрих мають по одному комплекту генетичної інформації; моногенні ознаки у таких організмів визначаються однією копією гена (на відміну від диплоїдних організмів, прикладом яких можуть бути люди, у яких навіть моногенні ознаки, як-от групи крові AB0, визначаються взаємодією двох копій відповідного гена; ці копії можуть бути однаковими, а можуть і відрізнятися);

• моноциклічні — такі організми, що розмножуються один раз протягом життя (і, звичайно, гинуть після цього); це не такий рідкісний варіант; наприклад, моноциклічними є усі головоногі молюски, до яких належать найбільші безхребетні та тварини з дуже складною поведінкою;

• клональні — тобто такі організми, розвиток яких проходить без рекомбінації, тобто без утворення нових поєднань спадкової інформації; при клональному відтворенні потомки генетично не відрізняються від батьківської особини;

• синхронне розмноження — не обов’язково таке, що відбувається одночасно, але таке, що можна моделювати для усіх особин, що розмножуються, на одному етапі циклу роботи моделі. Люди з цієї точки зору розмножуються асинхронно (їх розмноження може відбуватися на будь-якому етапі річного циклу та у особин різного віку), а, припустимо, пшениця — синхронно (навіть якщо на різних полях її розмноження може відбуватися у різні дні).

Зверніть увагу: ми не вважаємо синонімами поняття “відтворення” та “розмноження”. Відтворення відбувається на популяційному рівні біологічних систем. Розмноження відбувається на рівні організмів. Відтворення не прив’язано до певного часу року або етапу життєвого циклу; воно відбувається тому, що на заміну одним організмам приходять інші. Розмноження відбувається на певних етапах життєвого циклу та у певний час для кожного організму. Відтворюються популяції, а організми — розмножуються! Процеси організмового рівня, що забезпечують відтворення — це виживання та розмноження.

Тепер, після обговорення понять, ми можемо вказати, що за допомогою обговорюваної моделі ми хочемо встановити, якою буде динаміка двох форм організмів, що відрізняються за однією ознакою, у популяції з зазначеними властивостями за певних:

• початкових чисельностей обох форм;

• тією або іншою чисельністю популяції;

• під дією добору (певною перевагою однієї форми над іншою).

Сама по собі проблема, яку ми збираємося моделювати, є добре дослідженою за допомогою аналітичного моделювання. Існують формули, що дозволяють розрахувати очікувану чисельність на будь-який момент часу у разі відомих початкових параметрів. Втім, не завжди чисельність форм відповідатиме очікуваній. Чим менша популяція, тим більшим є внесок у кінцевий результат випадковості (“везіння” або “невдачі” однієї з форм). Якщо нас цікавить не найімовірніше очікуване значення, а, припустимо, розподіл ймовірностей різних значень, імітаційне моделювання дозволить отримати відповідь набагато зрозумілішим способом.

Є й інша перевага імітаційного моделювання. Імітаційну модель набагато простіше ускладнювати. Ми можемо додавати до обраної нами задачі додаткові умови, додавати відповідні розрахунки у модель та переходити від тривіальних до недостатньо вивчених випадків. Така перспектива можливого ускладнення надає сенсу рішенню поставленої задачі (і, навпаки, якби ми не планували ускладнювати обговорювану модель, її не було б сенсу будувати).

11.1.2 Відстежувати кожну особину окремо або просто підраховувати кількість однакових особин?

Це — дуже важливе питання. Конкретизуємо його.

Перший шлях: імітування індивідуальних “доль”. Кожній особині відповідає певне значення, яке зберігає модель. Якщо особина належить до категорії, де певна частина загине або буде відібрана для розмноження, ми маємо вирішувати це для неї “персонально”. Припустимо, половина особин загине. Тоді для кожної особини має бути розраховане випадкове число, за значенням якого має бути вирішено, вижила вони або загинула.

Другий шлях: підрахунки для однорідних груп. Ми підраховуємо кількість однакових особин і приймаємо рішення для усієї сукупності. Наприклад, якщо з певної кількості виживає половина, ми можемо розділити її на два (найпростіший варіант), або розрахувати кількість особин, що вижили, за розподілом Пуассона для даної кількості незалежних подій з заданою ймовірністю. Саме так були побудовані моделі експоненційного та логістичного зростання, які ми обговорювали у 4-му розділі.

Другий шлях має беззаперечну перевагу у швидкості виконання розрахунків. При достатньо великій чисельності особин розрахунки, що йдуть першим шляхом, надовго займуть навіть потужний комп’ютер, а за другим шляхом будуть виконані практично миттєво. Таким чином, використовувати другий шлях?

Не завжди. Сила імітаційного моделювання — у покроковому описі подій, що відбуваються з системою-оригіналом. Аналітичне моделювання майже завжди обирає другий шлях. Імітаційне використовує інший підхід; можливо, ключові переваги імітаційного моделювання більш наочно проявляться на першому шляху.

Врахуйте ще одну обставину. Ми починаємо з побудови дуже простої моделі. В ній усього дві групи особин. Чисельність популяції бути значною, і нам доведеться проводити розрахунки для N об’єктів; в розрахунках за другим шляхом ми матимемо справу з двома групами. Виграш другого шляху з точки зору кількості розрахунків беззаперечний. Але під час ускладнення моделі ми можемо дійти до того, що кожна або майже кожна особина в модельній популяції буде унікальною! Це знівелює переваги другого шляху розрахунків.

Студенти можуть самостійно використати других шлях розрахунків та порівняти результати з тими, що будуть отримані на першому шляху.

11.1.3 Послідовність розрахункових величин у робочому циклі моделі

З тієї задачі, яку ми поставили, зрозуміло, що ми моделюємо популяцію P, що складається з N особин. Ці особини не мають статі (саме так; вивчення еволюції статі ми починаємо з особин, що її позбавлені). Особини можуть належати до однієї з двох форм, які ми позначимо A та a. Ці форми — це різні генотипи (тобто варіанти спадкової інформації), яким відповідають різні фенотипи (результати розвитку особини, на які впливає спадкова інформація). Те, що ці дві форми відрізняються за фенотипами, ми знаємо на підставі того, що на їх співвідношення впливає добір. Добір завжди діє на фенотипи! Припустимо, добір сприяє формі a.

Популяція в цілому P на кожному етапі кожного робочого циклу імітаційної моделі (time) t є сумою двох субпопуляцій:
_tP^A + _tP^a = _tP.
Чисельність популяції (N) є сумою чисельностей двох форм:
_tN^A + _tN^a = _tN.
Внаслідок цього важливою характеристикою досліджуваної популяції є співвідношення чисельностей двох форм. Використаємо для його оцінки відносну чисельність форми a:
_tC^a = _tN^a / _tN
(позначення пов’язано з тим, що відносну чисельність форми у популяції можна розглядати як її концентрацію, concentration). Загалом позначення відповідають системі, яку ми обговорили у 4-мо розділі.

Нам треба відбити у моделі розмноження та добір. Розмноження передбачає певний життєвий цикл. Розглянемо найпростіший варіант: на кожному кроці моделі особини розмножуються і гинуть; до наступного кроку вони досягають зрілості і будуть здатні до розмноження знову.

Крім розмноження ми мали відбити у моделі добір. Що має відбуватися раніше: розмноження або добір? Для організмів, що розмножуються один раз, розмноження — результат їх життєвого циклу; події, пов’язані з добором, відбуваються раніше. Втім, цій логіці відповідає два варіанти організації циклу моделі.

Перший: початковий склад на початок циклу → розмноження → добір (загибель деяких особин) → остаточний склад на кінець циклу.

Другий: початковий склад на початок циклу → добір (загибель деяких особин) → розмноження → остаточний склад на кінець циклу.

Який з двох варіантів кращий? Перший. Він більш наочний: події, що стосуються однієї генерації (покоління) відбуваються протягом одного циклу роботи моделі. Можна врахувати й те, що у цьому випадку перший цикл можна почати з початкового складу особин (його слід задавати в початкових умовах), і саме особини цього початкового складу будуть розмножуватися у першому циклі роботи моделі.

Таким чином ми отримуємо низку зі станів модельної популяції, що послідовно змінюють один одного у ході роботи. Позначимо їх за допомогою грецьких літер.

• ^αP — альфа-сукупність — початковий склад плідників на кожному циклі роботи моделі; його чисельність — ^αN;

• ^βP — бета-сукупність — склад потомків (після розмноження); його чисельність — ^βN = ^αN × r;

• ^γP — гамма-сукупність — результат добору серед потомків;

• ^ωP — омега-сукупність — остаточний склад особин на кінець циклу.

З наведеного переліку етапів незрозуміло, чим гамма-сукупність відрізняється від омега-сукупності. Тут відбилися особливості реалізації моделі в R, які ми пояснимо дещо пізніше. Обговорювану модель можна реалізувати і так, щоб омега-сукупність розраховувалася відразу на підставі бета-сукупності.

11.1.4 Імітування розмноження

Як моделювати розмноження? Ми можемо використати найпростішій варіант. Згідно з експоненціальною моделлю зростання чисельності популяції, яку ми вже детально розглядали, _t+1N = r × _tN,
де r — репродуктивний потенціал, або ж параметр Мальтуса. Він характеризує здатність чисельності популяції до зростання у умовах, коли її ніщо не обмежує. Розрахунки у нашій моделі мають бути більш складними, оскільки у задачі моделювання входить аналіз добору, тобто певного обмеження зростання чисельності популяції.

Наведене рівняння стосується кількості особин (яку ми позначаємо N). Втім, аналогічну формулу можна написати й для популяції, як сукупності особин (P). Напишемо її не для наступного циклу (t+1), а для наступного етапу (бета-сукупності, а не альфа-сукупності). Таким чином,
^β_tP = r × ^α_tP.
Чи має такий вираз чіткій математичний сенс? Скоріше за все, ні, але він відповідає способу запису операцій з векторами у мові R (ви ж пам’ятаєте, що R — векторизована мова?).

Чи можуть репродуктивні потенціали двох форм організмів бути різними? Можуть. Якщо ми моделюємо розмноження та добір незалежно (це — не єдине можливе рішення, це можна було робити й одночасно), ми залишимо значення r однаковим для обох порівнюваних форм.

11.1.5 Функція sample і її використання для моделювання добору

Перед тим, як обговорити наступні два етапи роботи моделі, слід добре розібратися у роботі функції R, на використанні якої буде побудована реалізація цих етапів. Мова йде про функцію sample.

Функція sample(x, size, replace = FALSE, prob = NULL) може використовуватися як для генерації наборів випадкових цілих чисел, так і для вибору певної кількості об’єктів з їх сукупності. Два перші аргументи є обов’язковими: x — це діапазон, з якого здійснюється вибір, а size — розмір вибірки. Два наступних параметри є опціональними. Якщо вони не вказані, R використовує ті значення, що наведені у наведеному зразку повного запису аргументів. Так, параметр replace визначає, чи можна повторно обирати елемент, який вже був вибраний. У разі replace = FALSE після того, як обрано перший елемент, обране значення забирається з числа тих, які можуть бути обраними. Якщо вказати replace=TRUE, команда обиратиме кожний наступний елемент з повної сукупності.

Цікавим є параметр prob = NULL, який має задавати вектор з вагами розподілу (тобто задавати різну ймовірність вибору для різних елементів). За допомогою цього параметру можна було б імітувати і дію добору; втім, ми не будемо використовувати цю можливість, оскільки це зробить роботу нашої моделі менш наочною.

Розглянемо як працює функція sample на простому прикладі.

sample(1:10, 9)

## [1]  4 10  3  2  1  5  8  9  7

Вираз 1:10 задає вектор цілих чисел від одного до 10. З цієї сукупності береться 9 випадкових об’єктів. У результаті ми отримали сукупність з випадково розташованих 9 чисел з 10 (відсутнім лишилося якесь одно). Спробуйте виконати команду sample(1:10, 11)!

Чому R повідомляє про помилку? Якщо кожен об’єкт обирати один раз, неможливо обрати 11 з 10! Але якщо дозволити обирати числа повторно, команда буде виконана.

sample(1:10, 11, replace=T)

##  [1]  8  6 10  8  4  9  5  4  1  7  7

За допомогою цієї функції можна не лише генерувати випадкові числа, а й обирати об’єкти з певного переліку. У R є перелік букв англійської абетки, об’єкт letters. Створимо невипадкову вибірку з цього об’єкту, а потім — випадкову вибірку з цієї вибірки. Можна було обрати випадкові букви і з самого файлу letters.

dozen <- letters[1:12]
dozen

##  [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l"

sample(dozen, 4)

## [1] "i" "e" "j" "c"

sample(letters, 4)

## [1] "n" "b" "e" "z"

Ми можемо використати цю команду, щоб отримати випадкову вибірку, елементи в якій будуть з’являтися з різною ймовірністю. Припустимо, нам треба визначити “долю” 25 особин, кожна з яких має ймовірність вижити 70%, і, відповідно, ймовірність загинути 30%. Припустимо, ці особини перелічені у векторі Population.

Population <- 1:25
Population

##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

Кожна з цих цифр — це особина. Як визначити, якій жити, а якій помирати? Створимо вектор Fate, який задасть ймовірність вижити або померти.

Fate <- c(rep(NA, 3), rep(1, 7))
Fate

##  [1] NA NA NA  1  1  1  1  1  1  1

Символи, що визначають долю, розташовані у векторі Fate впорядковано, але якщо зробити з нього вибірку (з можливістю повторного обрання вже обраних елементів), цей порядок стане випадковим.

Lot <- sample(Fate, length(Population), replace=T)
Lot

##  [1]  1  1  1  1  1  1  1  1 NA NA  1  1  1  1  1  1  1  1  1  1  1 NA  1  1  1

Залишилося застосувати обраний жереб для визначення долі кожної особини. І тут ми використаємо особливості R як векторизованої мови. Річ у тім, що в R добуток векторів, якщо вони мають однакову довжину, це вектор тієї самої довжини, де перше значення — добуток перших значень векторів, друге — добуток других тощо.

Survivors <- Population*Lot
Survivors

##  [1]  1  2  3  4  5  6  7  8 NA NA 11 12 13 14 15 16 17 18 19 20 21 NA 23 24 25

Лишилося прибрати відсутні значення NA, що можна зробити, наприклад, командою na.omit(), і ми отримаємо скорочений склад популяції. До речі, цей перехід з прибиранням NA і відрізняє гама-чисельність від омега-чисельності і нашій моделі.

11.1.6 Імітування добору

Як моделювати добір? Для цього слід виміряти його інтенсивність. Звернемося до класичних визначень:

“Селективну перевагу одного алеля перед альтернативним алелем (або алелями) можна виразити у відсотках або у вигляді коефіцієнту добору (s), величина якого змінюється в діапазоні від 0 до 1. Кількісне значення коефіцієнта добору виводиться із відносних темпів відтворення альтернативних алелей. Припустимо, що в якійсь великій популяції а — алель, що має перевагу, а А — алель, якому добір не сприяє. У цій популяції на кожні 100 алелей а, переданих наступному поколінню, буде передаватися також кілька алелей А (від 100 до 0). Коефіцієнт добору є функцією цього відношення. Величину s можна визначити за формулою
s = 1 – (Швидкість відтворення алеля, якому не сприяє добір/Швидкість відтворення алеля, якому сприяє добір)”
Грант В. Эволюционный процесс, 2008

На успішність репродукції алелей впливає низка факторів. Носії різних алелей можуть відрізнятися за смертністю на різних етапах життєвого циклу, за різною імовірністю утворити пару у разі статевого розмноження (вона може відбиватися у різній привабливості для партнерів або різній ймовірності добитися їх у змаганні з конкурентами), різній плодючості тощо. Ми можемо штучно обмежити добір у нашої моделі одним моментом робочого циклу.

Як зімітувати добір з коефіцієнтом, припустимо, 0.3? У найпростішому для моделювання випадку — зберігати усіх особин, що належать до підтримуваної добором форми, і лише 70% (тобто 1–s) особин з форми, що не підтримується добором. Чому так?

Повернемося до формул, за допомогою яких ми дали визначення s:
Швидкість відтворення алеля, якому не сприяє добір = ^ωN^A/^αN^A
Швидкість відтворення алеля, якому сприяє добір = ^ωN^a/^αN^a

Ви пам’ятаєте, що
^αN × r = ^βN.
Якщо ми розділимо обидва наведені вирази на r, їх співвідношення не зміниться. У разі, якщо усі a-особини будуть переходити з бета-сукупності до гамма-сукупності,
^γN^a/^βN^a = 1.
У такому разі, якщо
^γN^A/^βN^A = 0.7, s = 1 – 0.7 = 0.3.

Ми використаємо змішану форму запису формул: вставимо у команди R не назви об’єктів, які будуть використані у справжній моделі, а позначення, які ми використовуємо для опису цієї моделі.

Припустимо, ми визначатимемо s з точністю до 0.01. Тоді ми можемо написати:
Fate <- c(rep(NA, s×100), rep(1, (1-s)×100))
Lot <- sample(Fate, length(^βN^A), replace=T)
^γN^A <- ^βN^A × Lot
^γN^a <- ^βN^a.

11.1.7 Імітування неселективного скорочення при виборі плідників для наступного циклу

При формуванні гамма ми використали добір, тобто диференційне виживання. Добір скорочує чисельність модельної популяції, а розмноження — збільшує. Ці два процеси не можуть бути точно врівноважені. Щоб модельна популяція була стійкою, розмноження має додавати більшу кількість особин, ніж віднімає добір. А як тоді запобігти розростанню популяції? За допомогою неселективного скорочення. Ми будемо проводити його при переході від омега-сукупності одного циклу до альфа-сукупності наступного. Нам треба вибрати з омега-чисельності попереднього циклу таку кількість особин для альфа-чисельності наступного циклу, яка буде дорівнювати N.

Як імітувати неселективне скорочення? За допомогою тієї ж самої функції sample(x, size, replace = FALSE, prob = NULL), яку можна використовувати без зайвих атрибутів: sample(x, size). Зрозуміло, що при формуванні сукупності плідників кожну особину з омега-сукупності можна обрати лише один раз:
^α_tP = sample(^ω_t–1P, N).

11.1.8 Алгоритм розрахунків

• ENTRANCE — ВХІД:

  • Initial state of the system — початковий стан модельованої системи

    • початкові чисельності обох форм: ^α₁N^A та ^α₁N^a;

  • Parameters — параметри

    • чисельність популяції: N;

    • параметр Мальтуса: r;

    • коефіцієнт добору s, де s ≈ 1 – (_t+1N^A/_tN^A) / (_t+1N^a/_tN^a).

  • Experimental conditions — умови експерименту з моделлю

    • кількість циклів моделі: cycles;

    • (за необхідності, у разі моделі ІІ типу): кількість ітерацій (прогонів моделі з однаковими початковими умовами) iterat;

  • The changeable parameters combinations — комбінації змінюваних параметрів

    • (за необхідності, у разі моделі ІІІ типу): закономірності перебору різних сполучень початкових параметрів.

• CALCULATIONS — РОЗРАХУНКИ:

  • Альфа-сукупність (плідники):

    • Якщо t = 1, чисельність субпопуляцій задається початковими параметрами, ^α₁N^A та ^α₁N^a;

    • Якщо t > 1, альфа-сукупність формується у ході неселективного скорочення омега-сукупності: ^α_tP = sample(^ω_t–1P, N).

  • Бета-сукупність (потомство):

    • ^β_tP = r × ^α_tP.

  • Гамма-сукупність (потомство, що проріджене добором):

    • ^γN^a <- ^βN^a
      Fate <- c(rep(NA, s×100), rep(1, (1-s)×100))
      Lot <- sample(Fate, length(^βN^A), replace=T)
      ^γN^A <- ^βN^A × Lot.

  • Омега-сукупність (з якої обираються плідники):

    • ^ωP = na.omit(^βP).

Насправді, модель вже створена; усе інше — реалізація викладеного тут задуму.

11.2.1 Початкові параметри та умови експерименту

Відкриваємо RStudio, створюємо новий файл (скрипт R) і починаємо писати модель! Спочатку задаємо інформаційний від у модель: початковий стан, параметри та умови експерименту з моделлю. Відкіля брати значення, які вказувати у моделі? Відкіля завгодно, “зі стелі”! Щоб писати скрипт та спостерігати його роботу, потрібно задати якісь значення. Наскільки вони змістовні, як вони співвідносяться з ціллю модельного експерименту, вирішувати треба буде пізніше. Спочатку краще обирати невеликі величини, щоб простіше було побачити у консолі проміжні результати обчислень.

# ДОБІР ЗА ОДНІЄЮ ОЗНАКОЮ В ПОПУЛЯЦІЇ МОНОЦИКЛІЧНИХ КЛОНАЛЬНИХ ОРГАНІЗМІВ З СИНХРОННИМ РОЗМНОЖЕННЯМ (Модель І типу)
# 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 # Кількість робочих циклів моделі 

Задавати початкові параметри можна як конкретними числами, так і певними залежностями (формулами), як ми це зробили для ^α₁N^A. Якщо ми будемо імітувати добір на користь форми a, логічно почати з випадку, коли частка a у модельній популяції (її “концентрація”) є мінімальною. Цей рядок (де використовується параметр N) має бути розташований після того, як параметр N буде визначений.

Під час створення моделі бажано задавати відносно невеликі значення параметрів та кількостей циклів. Під час покрокового дописування моделі (як буде показано надалі) буває корисно перевіряти, що виходить на кожному кроці. У такому разі зручно, щоб об’єкти, які будуть виводитися у консоль, не були занадто великими.

11.2.2 Створення об’єктів, з якими буде працювати модель

Планування об’єктів, які будуть створені у моделі, слід починати “з кінця”. У якому вигляді буде здійснюватися вивід результатів роботи моделі? Як треба організувати дані?

Скоріше за все, головним результатом певної ітерації є динаміка “концентрації” однієї з форм у часі, тобто зміни _tC^a залежно від t. У разі однієї ітерації цю динаміку може зберігати один вектор, у разі низки послідовних ітерацій — матриця. Почнемо з того варіанту, який будемо ускладнювати далі (ми його обрали, коли задали у початкових умовах кількість ітерацій, більшу за одиницю).

# II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
Rezult <- rep(NA, cycles+1) # Створення вектора для запису підсумків
Rezult

## [1] NA NA NA NA NA NA NA

Перед тим, як починати робочий цикл моделі, бажано створити ті об’єкти, що будуть у ньому використовуватися. R загалом є мовою, що не дуже пристосована до роботи з циклами. Найповільніше ця мова обробляє цикли, де створюються нові об’єкти. Тому до того, як запускати робочий цикл моделі, бажано створити об’єкти, з якими вона буде працювати. Перезаписувати значення у готовому об’єкті — простіша задача, ніж створювати новий об’єкт. Втім, для одного об’єкту, для омега-складу популяції, ми зробимо виключення і будемо створювати його всередині циклу. Річ у тім, що його довжина для різних циклів буде різною.

Al_Pop <- rep(NA, N) # Створення файлу для початкового складу
Be_Pop <- rep(NA, N*r) # Створення файлу для результатів розмноження
Ga_Pop <- rep(NA, N*r) # Створення файлу для результатів дії добору

Лишилося створити ще один об’єкт: зразок, випадкова вибірка з якого визначатиме ймовірність загибелі та виживання при доборі A-особин.

Fate <- c(rep(NA, s*100), rep(1, (1-s)*100)) # Вектор, вибірка з якого буде визначати долю у доборі 
Fate

##   [1] NA NA NA NA NA NA NA NA NA NA  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
##  [26]  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
##  [51]  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
##  [76]  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1

Альфа-склад популяції на першому циклі задається інакше, ніж на усіх наступних. Тож можна винести його перше визначення попереду, до початку робочого циклу. Раз ця величина розраховується окремо, то і її запис у файл з результатами слід здійснити окремо. Як ви бачите, форму a ми позначили 1, а форму A — 0. Використаний спосіб працює лише у тому випадку, коли ₁N^a ≥ 1.

# III. CALCULATIONS — РОЗРАХУНКИ:
# Initial composition creation — утворення початкового складу
Al_Pop[1:al_1_N_a] <- 1 # Перенесення до початкового складу a-особин
Al_Pop[(al_1_N_a+1):(al_1_N_a+al_1_N_A)] <- 0 # Перенесення до початкового складу A-особин
Al_Pop

##  [1] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Чи завжди треба виводити у консоль створені об’єкти? Звісно, ні. Але саме зараз, коли ми крок за кроком створюємо скрипт, це дуже корисно: можна відразу подивитися, чи отримуємо ми те, що треба. Під час написання скрипту дуже легко зробити якісь помилки; постійний зворотній зв’язок допомагає їх швидко знаходити та виправляти.

Якщо альфа-сукупність першого циклу створено, можна записати її характеристику у файл з результатами

Rezult[1] <- length(which(Al_Pop==1)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==0))) # Запис до вектора з підсумками першої позиції: початкового складу на початку роботи моделі
Rezult

## [1] 0.05   NA   NA   NA   NA   NA   NA

Як розраховувати концентрацію певної форми? Ми використали розрахунок “в лоб”, підрахували кількість a- та A-особин. У випадку позначень, що були використані, C^a можна було розрахувати й інакше, наприклад, так:

al_C_a <- sum(Al_Pop)/length(Al_Pop)
al_C_a

## [1] 0.05

Як ви бачите, результат тій самий. Так чому ж ми використали більш складний варіант? Заради універсальності. Якщо ми будемо відстежувати, припустимо, динаміку трьох форм, стане неможливим визначити кількість кожної з них шляхом обчислення суми усіх значень. А тій спосіб, який використали ми, буде працювати у любому випадку.

11.2.3 Компоненти робочого циклу моделі

Об’єкти створені, можна починати робочий цикл. Поки ми його створюємо, нам бажано отримувати зворотній зв’язок: перевіряти, чи отримуємо ми те, що планували. Як це робити, поки цикл не дописаний до кінця? Запланувати місце, де треба буде починати цикл, але тимчасового його закрити “ґраткою”(#). Цикл має перебирати значення змінної t. У такому разі, позначимо її вручну! Коли допишемо цикл до кінця і він запрацює, приберемо ці тимчасові конструкції.

# Main work cycle — основний робочий цикл
# for (t in 1:cycles) { # Робочий цикл моделі
t <- 1 # після запуску циклу цю строку можна буде заблокувати або прибрати
# (тут будемо покроково створювати цикл)
#                     } # Кінець робочого циклу

Починаємо. Створимо вектор для потомства і подивимося, що вийшло.

Be_Pop <- rep(Al_Pop, r) # Заповнення вектора для потомства
Be_Pop

##  [1] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [39] 0 0

Те, що треба, так? Тепер нам треба визначити кількості особин обох форм (їх треба буде враховувати при заповненні вектора з гамма-чисельністю).

Be_N_a <- length(which(Be_Pop==1)) # Розрахунок кількості a-особин
Be_N_A <- length(which(Be_Pop==0)) # Розрахунок кількості A-особин
Be_N_a

## [1] 2

Be_N_A

## [1] 38

Можна починати заповнювати вектор Ga_Pop, що містить результати добору. Як обговорювалося раніше, a-особини безперешкодно потрапляють до гамма-сукупності.

Ga_Pop[1:Be_N_a] <- subset(Be_Pop, Be_Pop==1) # a-особини зберігаються усі

Жереб, який отримають A-особини, визначає вектор Lot.

Lot <- sample(Fate, Be_N_A, replace=T) # Вибірка, що визначає долю A-особин
Ga_Pop[(Be_N_a+1):(Be_N_a+Be_N_A)] <- subset(Be_Pop, Be_Pop==0)*Lot # Загиблі позначені NA
Ga_Pop

##  [1]  1  1  0 NA  0 NA  0 NA  0  0  0  0  0  0  0  0  0  0  0  0  0  0 NA  0  0
## [26]  0  0  0 NA  0  0  0  0  0  0  0  0 NA  0  0

Тепер можна отримати омега-чисельність, прибрав зі складу модельної популяції NA, які маркували загиблих особин.

Om_Pop <- na.omit(Ga_Pop) # Ті, хто вижили (без врахування загиблих)
Om_Pop

##  [1] 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## attr(,"na.action")
## [1]  4  6  8 23 29 38
## attr(,"class")
## [1] "omit"

Лишилося провести неселективне скорочення чисельності модельної популяції до N.

Al_Pop <- sample(Om_Pop, N) # Плідники наступного покоління
Al_Pop

##  [1] 0 0 0 0 0 0 0 1 0 0

Усі компоненти робочого циклу моделі задані? Ні. На попередньому кроці циклу t ми отримали альфа-склад модельної популяції на циклі t+1. Його слід записати у матрицю Rezult. Але перед тим слід додати ще два рядки, які спрацюють, якщо у модельній популяції залишиться лише одна форма.

Якщо в альфа-складі попереднього циклу були a-особини, вони мали потрапити до омега-складу. Але якщо їх небагато, вони у силу простої імовірності можуть не потрапити до альфа-чисельності наступного циклу. Чи є при цьому сенс проводити наступні розрахунки у даній ітерації? Ні, якщо в модельній популяції лишилися особини тільки однієї форми, склад цієї популяції змінюватися вже не буде. Ітерації. слід припинити і записати, що вона закінчилася зникненням a-особин. Де записати? У останньому рядку матриці Rezult!

Аналогічно слід передбачити і можливість зникнення з модельної популяції A-особин; вони могли зникнути і під час переходу від бета-складу до гамма складу, і при переході від омега-складу до альфа-складу. Щоб передбачити ці можливості, використаємо умови if (умова) {…; …}. У круглих дужках ставиться умова, що перевіряється. У разі, якщо ця умова відповідає дійсності, виконуються подальші команди. Низку команд, що мають бути виконані послідовно, можна об’єднати у фігурних дужках і або записати на послідовних рядках, або розділити крапками з комами. Функція break перериває виконання циклу for, тобто забезпечує переривання перебору послідовних робочих циклів і перехід до наступної ітерації.

При запуску наступного фрагменту R-коду R Markdown “не розуміє”, до якого циклу відносяться команди break, адже робочий цикл моделі зараз не виконується (ми його ще не ввімкнули). Закриємо рядки з цими командами ґратками.

# "Запобіжники" на випадок зникнення однієї з форм:
#  if (length(which(Al_Pop==1))==0) {Rezult[cycles+1] <- 0; break} # У разі перемоги A-особин 
#  if (length(which(Al_Pop==0))==0) {Rezult[cycles+1] <- 1; break} # У разі перемоги a-особин
# II.4 Збір підсумків:
  Rezult[t+1] <- length(which(Al_Pop==1)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==0))) # Запис наступного рядка до вектора з результатами роботи
  Rezult

## [1] 0.05 0.10   NA   NA   NA   NA   NA

#                     } # Кінець робочого циклу

Хоча робочий цикл моделі ще не запущений, у векторі Rezult на першій та другій позиції мали з’явитися відповідні записи.

11.2.4 Запуск робочого циклу моделі

Ви пам’ятаєте, що ми у нашому скрипті позначили місце початку робочого циклу? Тепер його слід запустити. Наведемо усю модель, в яку поки що не включено візуалізацію результатів. Використаємо такі вхідні умови, за якими буде цікавіше спостерігати за моделлю

# ДОБІР ЗА ОДНІЄЮ ОЗНАКОЮ В ПОПУЛЯЦІЇ МОНОЦИКЛІЧНИХ КЛОНАЛЬНИХ ОРГАНІЗМІВ З СИНХРОННИМ РОЗМНОЖЕННЯМ (Модель І типу)
# I. ENTRANCE — ВХІД:
# Initial script commands — початкові команди скрипту
# setwd("~/!_Courses/Sex_on_R") # Робоча директорія (лише на комп'ютері Д.Ш.!!!)
rm(list = ls()) # Очищення раніше збережених об'єктів в Environment
# Initial state of the system — початковий стан модельованої системи
al_1_N_a <- 1 # Початкова чисельність a-особин 
# Parameters — параметри
N <- 50 # Альфа-чисельність модельної популяції
al_1_N_A <- N - al_1_N_a # Початкова чисельність A-особин (усі крім a-особин)
r <- 2 # Параметр Мальтуса (кількість потомків на одного плідника на один цикл роботи моделі)
s <- 0.2 # Коефіцієнт добору на користь a-особин 
# Experimental conditions — умови експерименту з моделлю
cycles <- 50 # Кількість робочих циклів моделі
# II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
Rezult <- rep(NA, cycles+1) # Створення вектора для запису підсумків
Al_Pop <- rep(NA, N) # Створення файлу для початкового складу
Be_Pop <- rep(NA, N*r) # Створення файлу для результатів розмноження
Ga_Pop <- rep(NA, N*r) # Створення файлу для результатів дії добору
Fate <- c(rep(NA, s*100), rep(1, (1-s)*100)) # Вектор, вибірка з якого буде визначати долю у доборі 
# III. CALCULATIONS — РОЗРАХУНКИ:
# Initial composition creation — утворення початкового складу
Al_Pop[1:al_1_N_a] <- 1 # Перенесення до початкового складу a-особин
Al_Pop[(al_1_N_a+1):(al_1_N_a+al_1_N_A)] <- 0 # Перенесення до початкового складу A-особин
Rezult[1] <- length(which(Al_Pop==1)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==0))) # Запис до вектора з підсумками першої позиції: початкового складу на початку роботи моделі
# Main work cycle — основний робочий цикл
for (t in 1:cycles) { # Робочий цикл моделі
  Be_Pop <- rep(Al_Pop, r) # Заповнення вектора для потомства
  Be_N_a <- length(which(Be_Pop==1)) # Розрахунок кількості a-особин
  Be_N_A <- length(which(Be_Pop==0)) # Розрахунок кількості A-особин
  Ga_Pop[1:Be_N_a] <- subset(Be_Pop, Be_Pop==1) # a-особини зберігаються усі
  Lot <- sample(Fate, Be_N_A, replace = T) # Вибірка, що визначає долю A-особин
  Ga_Pop[(Be_N_a+1):(Be_N_a+Be_N_A)] <- subset(Be_Pop, Be_Pop==0)*Lot # Загиблі позначені NA
  Om_Pop <- na.omit(Ga_Pop) # Ті, хто вижили (без врахування загиблих)
  Al_Pop <- sample(Om_Pop, N) # Плідники наступного покоління
# "Запобіжники" на випадок зникнення однієї з форм:
  if (length(which(Al_Pop==1))==0) {Rezult[cycles+1] <- 0; break} # У разі перемоги A-особин 
  if (length(which(Al_Pop==0))==0) {Rezult[cycles+1] <- 1; break} # У разі перемоги a-особин
  Rezult[t+1] <- length(which(Al_Pop==1)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==0))) # Запис наступного рядка до вектора з результатами роботи
                     } # Кінець робочого циклу
# IV. FINISHING — ЗАВЕРШЕННЯ:
# Results viewing — перегляд підсумків
Rezult # Перегляд вектора з динамікою концентрації a-особин

##  [1] 0.02 0.02 0.04 0.06 0.06 0.10 0.12 0.16 0.26 0.26 0.34 0.50 0.56 0.62 0.70
## [16] 0.78 0.78 0.76 0.82 0.88 0.92 0.92 0.90 0.90 0.88 0.94 0.96 0.98   NA   NA
## [31]   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
## [46]   NA   NA   NA   NA   NA 1.00

Працює!!!

Лишилося додати графік

# Results visualization — візуалізація підсумків
plot(Rezult, type="l", lty=1, col="red") # Ця команда виводить результати на найпростішій графік

Модель (в її мінімальній комплектації) запрацювала. Але це ще далеко не усе.

11.2.5 Приклад роботи моделі

Запустимо програму у показаному далі варіанті (збільшимо чисельність популяції і кількість циклів, щоб подивитися на масштабнішу картину, ніж під час отладки). Зверніть увагу: команда set.seed() є активною, тобто раз за разом виконання програми буде приводити до того ж самого результату. Останній рядок — команда відбити на графіку результати першої ітерації. Виводити на консоль матрицю Rezult не будемо, але ви, якщо цікаво, можете подивитися її самі.

# ДОБІР ЗА ОДНІЄЮ ОЗНАКОЮ В ПОПУЛЯЦІЇ МОНОЦИКЛІЧНИХ КЛОНАЛЬНИХ ОРГАНІЗМІВ З СИНХРОННИМ РОЗМНОЖЕННЯМ (Модель І типу)
# I. ENTRANCE — ВХІД:
# Initial script commands — початкові команди скрипту
# setwd("~/!_Courses/Sex_on_R") # Робоча директорія (лише на комп'ютері Д.Ш.!!!)
rm(list = ls()) # Очищення раніше збережених об'єктів в Environment
set.seed(1234567)
# 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.01 # Коефіцієнт добору на користь a-особин 
# Experimental conditions — умови експерименту з моделлю
cycles <- 100 # Кількість робочих циклів моделі
# II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
Rezult <- rep(NA, cycles+1) # Створення вектора для запису підсумків
Al_Pop <- rep(NA, N) # Створення файлу для початкового складу
Be_Pop <- rep(NA, N*r) # Створення файлу для результатів розмноження
Ga_Pop <- rep(NA, N*r) # Створення файлу для результатів дії добору
Fate <- c(rep(NA, s*100), rep(1, (1-s)*100)) # Вектор, вибірка з якого буде визначати долю у доборі 
# III. CALCULATIONS — РОЗРАХУНКИ:
# Initial composition creation — утворення початкового складу
Al_Pop[1:al_1_N_a] <- 1 # Перенесення до початкового складу a-особин
Al_Pop[(al_1_N_a+1):(al_1_N_a+al_1_N_A)] <- 0 # Перенесення до початкового складу A-особин
Rezult[1] <- length(which(Al_Pop==1)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==0))) # Запис до вектора з підсумками першої позиції: початкового складу на початку роботи моделі
# Main work cycle — основний робочий цикл
for (t in 1:cycles) { # Робочий цикл моделі
  Be_Pop <- rep(Al_Pop, r) # Заповнення вектора для потомства
  Be_N_a <- length(which(Be_Pop==1)) # Розрахунок кількості a-особин
  Be_N_A <- length(which(Be_Pop==0)) # Розрахунок кількості A-особин
  Ga_Pop[1:Be_N_a] <- subset(Be_Pop, Be_Pop==1) # a-особини зберігаються усі
  Lot <- sample(Fate, Be_N_A, replace = T) # Вибірка, що визначає долю A-особин
  Ga_Pop[(Be_N_a+1):(Be_N_a+Be_N_A)] <- subset(Be_Pop, Be_Pop==0)*Lot # Загиблі позначені NA
  Om_Pop <- na.omit(Ga_Pop) # Ті, хто вижили (без врахування загиблих)
  Al_Pop <- sample(Om_Pop, N) # Плідники наступного покоління
# "Запобіжники" на випадок зникнення однієї з форм:
  if (length(which(Al_Pop==1))==0) {Rezult[cycles+1] <- 0; break} # У разі перемоги A-особин 
  if (length(which(Al_Pop==0))==0) {Rezult[cycles+1] <- 1; break} # У разі перемоги a-особин
  Rezult[t+1] <- length(which(Al_Pop==1)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==0))) # Запис наступного рядка до вектора з результатами роботи
                     } # Кінець робочого циклу
# IV. FINISHING — ЗАВЕРШЕННЯ:
# Results viewing — перегляд підсумків
Rezult # Перегляд вектора з динамікою концентрації a-особин

##   [1] 0.01 0.02 0.03 0.02 0.01 0.02 0.03 0.02 0.02 0.04 0.05 0.06 0.08 0.09 0.06
##  [16] 0.06 0.05 0.08 0.12 0.13 0.10 0.09 0.07 0.09 0.10 0.11 0.14 0.13 0.21 0.19
##  [31] 0.20 0.17 0.12 0.11 0.08 0.09 0.13 0.16 0.18 0.17 0.18 0.17 0.19 0.17 0.19
##  [46] 0.18 0.15 0.17 0.23 0.20 0.14 0.15 0.16 0.18 0.13 0.14 0.15 0.17 0.21 0.24
##  [61] 0.26 0.20 0.18 0.12 0.12 0.10 0.10 0.09 0.08 0.09 0.11 0.09 0.07 0.05 0.05
##  [76] 0.05 0.03 0.04 0.04 0.04 0.05 0.04 0.04 0.05 0.08 0.10 0.11 0.08 0.06 0.05
##  [91] 0.04 0.03 0.04 0.04 0.05 0.04 0.05 0.04 0.02 0.01 0.01

# Results visualization — візуалізація підсумків
plot(Rezult, type="l", lty=1, col="red") # Ця команда виводить результати на найпростішій графік

Що ми маємо? Частка (“концентрація”) a-особин, показана червоною лінією, протягом певного часу скоріше зростала, чим знижувалася, а потім скоріше знижувалася, ніж зростала. Скінчилося це тим, що усі a-особини зникли з модельної популяції. Таку динаміку досліджуваної величини можна назвати блуканням: вона визначається переважно випадковими обставинами, а не дією добору.

Так що, у таких умовах, які були використані у моделі, добір не може змінити склад модельної популяції? Може, але не завжди…

11.3.1 Перебудова моделі І типу у модель ІІ типу

Як ми зазначили у 6 розділі, моделі І типу проводять одну ітерацію (один прохід моделі з певними умовами), моделі ІІ типу проводять групу ітерацій з однаковими умовами, а моделі ІІІ типу проводять групи ітерацій, що відповідають різним значенням досліджуваних параметрів. У підрозділі 11.3 ми створили модель І типу. Що потрібно, щоб перетворити її на модель ІІ типу?

Необхідно додати засоби, що забезпечують проведення групи ітерацій з однаковими умовами. Це нескладно. Ми “натягнемо” цикл ітерацій на робочий цикл моделі.

Які зміни в моделі І типу для цього потрібні? Перш за все, якщо ми плануємо отримати об’єкт, що буде містити усі підсумки (зміну “концентрації” особин, яким сприяє добір), нам незручно буде мати справу з вектором. Зручніше зробити матрицю, де кожний стовпець відповідає певній ітерації.

У разі, якщо об’єкт Rezult перетворився на матрицю, посилання на нього мають бути двовимірними. Нам недостатньо вказати, до якого циклу (t) належить те значення, з яким ми маємо справу; нам треба також вказати, до якої ітерації належить це значення: Rezult[t, i].

В усьому іншому наша модель ІІ типу не відрізняється від моделі І типу.

# ДОБІР ЗА ОДНІЄЮ ОЗНАКОЮ В ПОПУЛЯЦІЇ МОНОЦИКЛІЧНИХ КЛОНАЛЬНИХ ОРГАНІЗМІВ З СИНХРОННИМ РОЗМНОЖЕННЯМ (Модель ІІ типу)
# I. ENTRANCE — ВХІД:
# Initial script commands — початкові команди скрипту
# setwd("~/!_Courses/Sex_on_R") # Робоча директорія (лише на комп'ютері Д.Ш.!!!)
rm(list = ls()) # Очищення раніше збережених об'єктів в Environment
set.seed(1234567)
# 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.01 # Коефіцієнт добору на користь a-особин 
# Experimental conditions — умови експерименту з моделлю
cycles <- 10 # Кількість робочих циклів моделі 
iterat <- 3 # Кількість ітерацій
# II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
Rezult <- matrix(NA, nrow = cycles+1, ncol = iterat) # У разі використання ітерацій, потрібна матриця, а не вектор, як у моделі І типу
# III. CALCULATIONS — РОЗРАХУНКИ:
# Higher level cycles running — запуск циклів вищого рівня
for (i in 1:iterat) { # Початок ітерації
  Al_Pop <- rep(NA, N) # Створення файлу для початкового складу
  Be_Pop <- rep(NA, N*r) # Створення файлу для результатів розмноження
  Ga_Pop <- rep(NA, N*r) # Створення файлу для результатів дії добору
  Fate <- c(rep(NA, s*100), rep(1, (1-s)*100)) # Вектор, вибірка з якого буде визначати долю у доборі 
  # Initial composition creation — утворення початкового складу
  Al_Pop[1:al_1_N_a] <- 1 # Перенесення до початкового складу a-особин
  Al_Pop[(al_1_N_a+1):(al_1_N_a+al_1_N_A)] <- 0 # Перенесення до початкового складу A-особин
  Rezult[1, i] <- length(which(Al_Pop==1)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==0))) # Запис до матриці з підсумками першої позиції: початкового складу на початку роботи моделі
  # Main work cycle — основний робочий цикл
  for (t in 1:cycles) { # Робочий цикл моделі
    Be_Pop <- rep(Al_Pop, r) # Заповнення вектора для потомства
    Be_N_a <- length(which(Be_Pop==1)) # Розрахунок кількості a-особин
    Be_N_A <- length(which(Be_Pop==0)) # Розрахунок кількості A-особин
    Ga_Pop[1:Be_N_a] <- subset(Be_Pop, Be_Pop==1) # a-особини зберігаються усі
    Lot <- sample(Fate, Be_N_A, replace = T) # Вибірка, що визначає долю A-особин
    Ga_Pop[(Be_N_a+1):(Be_N_a+Be_N_A)] <- subset(Be_Pop, Be_Pop==0)*Lot # Загиблі позначені NA
    Om_Pop <- na.omit(Ga_Pop) # Ті, хто вижили (без врахування загиблих)
    Al_Pop <- sample(Om_Pop, N) # Плідники наступного покоління
    # "Запобіжники" на випадок зникнення однієї з форм:
    if (length(which(Al_Pop==1))==0) {Rezult[cycles+1, i] <- 0; break} # У разі перемоги A-особин 
    if (length(which(Al_Pop==0))==0) {Rezult[cycles+1, i] <- 1; break} # У разі перемоги a-особин
    # Ending of higher level cycles — завершення циклів вищого рівня
    Rezult[t+1, i] <- length(which(Al_Pop==1)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==0))) # Запис наступного рядка до вектора з результатами роботи
                       } # Кінець робочого циклу
                     } # Кінець ітерації
# IV. FINISHING — ЗАВЕРШЕННЯ:
# Results viewing — перегляд підсумків
Rezult

##       [,1] [,2] [,3]
##  [1,] 0.01 0.01 0.01
##  [2,] 0.02 0.02 0.02
##  [3,] 0.03 0.02 0.03
##  [4,] 0.02   NA 0.03
##  [5,] 0.01   NA 0.03
##  [6,] 0.02   NA 0.03
##  [7,] 0.03   NA 0.04
##  [8,] 0.02   NA 0.05
##  [9,] 0.02   NA 0.05
## [10,] 0.04   NA 0.05
## [11,] 0.05 0.00 0.05

plot(Rezult[, 1], type="l", lty=1, col="red") # Побудова графіку за результатами першої ітерації

Ми використали команду set.seed(1234567) і тепер впевнені, що при кожному виконання файлу R Markdown результат буде тим самим. З трьох ітерацій друга завершилася зникненням особин, яким сприяє добір. На графік ми вивели результати першої ітерації.

У одному випадку з трьох a-особини зникли… А якою є ймовірність їхньої перемоги? Для цього нам треба проаналізувати більшу за кількістю групу ітерацій, кожна з яких складатиметься з достатньої кількості циклів, щоб стало зрозуміло, чим вона закінчилася…

11.3.2 Аналіз групи ітерації за допомогою моделі ІІ типу

Запустимо модель ІІ типу з вказаними параметрами:

# ДОБІР ЗА ОДНІЄЮ ОЗНАКОЮ В ПОПУЛЯЦІЇ МОНОЦИКЛІЧНИХ КЛОНАЛЬНИХ ОРГАНІЗМІВ З СИНХРОННИМ РОЗМНОЖЕННЯМ (Модель ІІ типу)
# I. ENTRANCE — ВХІД:
# Initial script commands — початкові команди скрипту
# setwd("~/!_Courses/Sex_on_R") # Робоча директорія (лише на комп'ютері Д.Ш.!!!)
rm(list = ls()) # Очищення раніше збережених об'єктів в Environment
set.seed(1234567)
# 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.01 # Коефіцієнт добору на користь a-особин 
# Experimental conditions — умови експерименту з моделлю
cycles <- 300 # Кількість робочих циклів моделі 
iterat <- 100 # Кількість ітерацій
# II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
Rezult <- matrix(NA, nrow = cycles+1, ncol = iterat) # У разі використання ітерацій, потрібна матриця, а не вектор, як у моделі І типу
# III. CALCULATIONS — РОЗРАХУНКИ:
# Higher level cycles running — запуск циклів вищого рівня
for (i in 1:iterat) { # Початок ітерації
  Al_Pop <- rep(NA, N) # Створення файлу для початкового складу
  Be_Pop <- rep(NA, N*r) # Створення файлу для результатів розмноження
  Ga_Pop <- rep(NA, N*r) # Створення файлу для результатів дії добору
  Fate <- c(rep(NA, s*100), rep(1, (1-s)*100)) # Вектор, вибірка з якого буде визначати долю у доборі 
  # Initial composition creation — утворення початкового складу
  Al_Pop[1:al_1_N_a] <- 1 # Перенесення до початкового складу a-особин
  Al_Pop[(al_1_N_a+1):(al_1_N_a+al_1_N_A)] <- 0 # Перенесення до початкового складу A-особин
  Rezult[1, i] <- length(which(Al_Pop==1)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==0))) # Запис до матриці з підсумками першої позиції: початкового складу на початку роботи моделі
  # Main work cycle — основний робочий цикл
  for (t in 1:cycles) { # Робочий цикл моделі
    Be_Pop <- rep(Al_Pop, r) # Заповнення вектора для потомства
    Be_N_a <- length(which(Be_Pop==1)) # Розрахунок кількості a-особин
    Be_N_A <- length(which(Be_Pop==0)) # Розрахунок кількості A-особин
    Ga_Pop[1:Be_N_a] <- subset(Be_Pop, Be_Pop==1) # a-особини зберігаються усі
    Lot <- sample(Fate, Be_N_A, replace = T) # Вибірка, що визначає долю A-особин
    Ga_Pop[(Be_N_a+1):(Be_N_a+Be_N_A)] <- subset(Be_Pop, Be_Pop==0)*Lot # Загиблі позначені NA
    Om_Pop <- na.omit(Ga_Pop) # Ті, хто вижили (без врахування загиблих)
    Al_Pop <- sample(Om_Pop, N) # Плідники наступного покоління
    # "Запобіжники" на випадок зникнення однієї з форм:
    if (length(which(Al_Pop==1))==0) {Rezult[cycles+1, i] <- 0; break} # У разі перемоги A-особин 
    if (length(which(Al_Pop==0))==0) {Rezult[cycles+1, i] <- 1; break} # У разі перемоги a-особин
    # Ending of higher level cycles — завершення циклів вищого рівня
    Rezult[t+1, i] <- length(which(Al_Pop==1)) / (length(which(Al_Pop==1)) + length(which(Al_Pop==0))) # Запис наступного рядка до вектора з результатами роботи
                       } # Кінець робочого циклу
                     } # Кінець ітерації
# IV. FINISHING — ЗАВЕРШЕННЯ:
# Results viewing — перегляд підсумків
Rezult[cycles+1, ]

##   [1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
##  [16] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
##  [31] 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00
##  [46] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00
##  [61] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
##  [76] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.78 0.00
##  [91] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00

mean(Rezult[cycles+1, ]) # Перегляд матриці з динамікою концентрації a-особин протягом кожної ітерації

## [1] 0.0578

# plot(Rezult[, 1], type="l", lty=1, col="red") # Побудова графіку за результатами першої ітерації

Переконайтеся самі: у разі запуску моделі з вказаними нижче параметрами та умовами 5 ітерацій із 100 закінчуються перемогою a-особин, а в одній ітерації блукання ще триває (результат змагання a- та A-особин не визначений).

Може, s=0.01 — це такий слабкий добір, що він не впливає на модельну популяцію?

Це дійсно слабкий добір, але він усе ж таки впливає на популяцію. Що переконатися в цьому, задайте s <- 0 і запустіть модель з усіма іншими вказаними параметрами, що відповідають попередньому випадку. Ви побачите: a-особини перемогли лише у одному випадку, і ще в одній ітерації блукання тривало. Результат гірший, ніж з добором, хоча б слабким.

А тепер встановіть s <- 0.1: a-особини отримають 26 перемог зі 100 можливих! Втім, у 74 випадках вони усе-таки програють. Дія добору проявиться ще у тому, що результати ітерацій будуть визначатися швидше, і довге блукання стане нехарактерним.

Ми можемо дійти до висновку, що результатом роботи моделі ІІ типу є вже не результати окремої ітерації (вони можуть суттєво відрізнятися внаслідок випадкових причин), а розподіл ймовірностей результатів.

Як порівнювати розподіли ймовірностей результатів? Один з шляхів — так, як ми тільки-но робили. Змінювати параметри, вплив яких ми бажаємо встановити, та порівнювати результати різних груп ітерацій. Ми робили це вручну, і вручну характеризували отриману різницю між групами ітерацій. Але ж ми маємо справу з мовою програмування! Зробимо це із застосуванням її засобів.