?

Log in

No account? Create an account

Tue, Jan. 11th, 2011, 08:25 am
Сам я водопроводчик, но думаю, что роды надо принимать вот так…

Сам я водопроводчик, но думаю, что роды надо принимать вот так…

       Я - неисправимый фантазер. Только, в отличие от Поршнева, (http://victor-chapaev.livejournal.com/5439.html) стараюсь фантазировать в рамках реальных фактов. Больше скажу, без фактов и фантазировать неинтересно. На этот раз фантазию разбудила статья уважаемого Александра Маркова с коллегами под названием «Взаимосвязь размера генома и сложности организма в эволюционном ряду от прокариот к млекопитающим» в «Палеонтологическом журнале» № 4 за 2010 год, любезно размещенная им на своем сайте «Проблемы эволюции».

http://evolbiol.ru/large_files/mgenome2010.pdf

Коротко, о чем статья. Авторы изучали вопрос – есть ли корреляция между размером генома организмов и сложностью этих организмов. Теоретически рассуждая, такая корреляция должна быть – чем сложнее организм, тем больше нужно генетического материала, чтобы закодировать и передать потомству эту сложность. Но экспериментальные данные говорят о другом - наибольший размер геномов, например, обнаружен у относительно примитивных одноклеточных эукариотов (амёб), а многие земноводные имеют геномы, какие и не снились ни одному млекопитающему.

Подход авторов был таков: размер генома – компромисс между процессами увеличения генома и процессами его сокращения. Первые – это случайное удвоение части генетического материала из-за ошибок копирования, предположительное заимствование части чужого генетического материала в результате горизонтального переноса и т.п. Вторые – делеции части генома из-за тех же ошибок копирования. Причем и те и другие процессы, являясь, по сути, мутациями, могут быть смертельны, безразличны или полезны для организма. Авторы приняли во внимание, что геном любого, сколь угодно простого организма, может относительно безболезненно увеличится до сколь угодно большого размера, не нанося владельцу большого урона – увеличатся лишь расходы на его копирование. А вот уменьшится до сколь угодно малого размера геном не может – в конце концов, уменьшение затронет жизненно-важные части генома и мутация станет смертельной. Поэтому, если среди группы видов найти вид с наименьшим размером генома, он и будет характеризовать минимально-необходимое количество генетического материала, достаточное, чтобы закодировать сложность, характерную для этой группы.

Авторы выбрали девять вложенных друг в друга групп организмов: все живое, эукариоты (имеющие клеточные ядра и органеллы), животные, хордовые, позвоночные, тетраподы (четырехногие), амниоты (объединяют пресмыкающихся, птиц и млекопитающих), собственно млекопитающих и, наконец, приматов. Для каждой группы был определен примерный возраст ее появления на лице матушки-Земли и минимальный размер генома (в нуклеотидных парах). Фотография таблицы прилагается.


Рисунок 1.

 

Далее авторы от биологии переходят к математике, точнее, к матстатистике. Положив точки на график, авторы пытаются подобрать математическую функцию, или, другими словами, закон изменения размера генома от времени, который наилучшим образом мог бы описать весь набор данных. По оси ординат, по которой откладывается размер генома (в млн. пар нуклеотидов), выбран логарифмический масштаб, так как за 4 миллиарда лет, по данным авторов, минимальный размер генома вырос на три с половиной порядка.

Было испробовано два типа законов – двух- и трехпараметрические. А так как в логарифмическом масштабе наблюдается ускоренный рост (теперь уже даже не размера генома, а логарифма размера генома), то и выбирались функции нелинейные, дающие этот ускоренный рост. Из двухпараметрических моделей таковыми являлись рост по экспоненте (y=A*exp(-B*x)) и рост по гиперболе (y=A/(B+x)). Из трехпараметрических – степенной экспоненциальный рост (y=A*exp(-B*xN)), степенной гиперболический рост (y=A/(B+x)N) и двойной экспоненциальный (y=C*exp(A*exp(-B*x))). Коэффициенты A, B, C и N подбирались, как написано в статье, «прямым численным перебором при помощи метода наименьших квадратов». В результате проделанной вычислительной работы в каждом классе функций (двух- и трехпараметрических) выявлены наилучшие кандидаты. Ими оказались в классе двухпараметрических функций – закон экспоненциального роста, в классе трехпараметрических – двойной экспоненциальный рост. Приведены параметры этих аппроксимаций и значения «критерия Пирсона R2».

Данные доступны, попробую и я. Вот повторение результатов из статьи (рисунок 2). Сначала я, памятуя, что в авторской работе подбор параметров осуществлялся численным перебором, решил идти тем же путем. Только на домашнем ноутбуке, да и на рабочей станции на работе, простой перебор занял бы слишком много времени (особенно для трех параметров). Поэтому я, вспомнив, что в любимой мною программе MS Excel, входящей в стандартную поставку MS Office, имеется сервис под названием «поиск решения», применил именно его. Этот сервис как раз и производит поиск коэффициентов методом перебора, только не простого, а оптимизированного, близкого к методу градиентного спуска. Для трех параметров это миллисекундное дело. Завел столбцы исходных данных из статьи – x и y. Задал начальные значения коэффициентов для закона экспоненциального роста. Ввел расчетную формулу для каждого x. Посчитал разницы между логарифмом каждого игрека и логарифмом его аппроксимации по формуле, их квадрат и сумму этих квадратов. Задал поиск минимума суммы квадратов, область, где лежат коэффициенты, и запустил сервис. Полученный ответ (значения коэффициентов уравнения) совпал с приведенными в статье до долей процента. Посчитал критерий Пирсона (тоже стандартная функция Excel), построил график. График мой и построенный по коэффициентам из статьи визуально не различались (см. рис.2). Зато критерий Пирсона чего-то не совпал. В статье написано «R2= 0.914», а у меня КВПИРСОН()=0.8788, корень квадратный из него - 0.9375. Тогда посчитал КВПИРСОН() для логарифмов значений и аппроксимаций – вот оно: КВПИРСОН()=0.8364, а корень из него как раз и совпадает с приведенным в статье значением – 0.9145. Уф, видно, опечатка. Результаты усилий - на графике, там же показан график двойного экспоненциального роста по коэффициентам из статьи:


Рисунок 2.

 

В попытке получить коэффициенты двойного экспоненциального роста не преуспел, алгоритм где-то расходится. В зашкал уходят коэффициенты… А потом сообразил, что все эти формулы представляются простыми зависимостями - если построить график не самих значений размера генома, а их логарифмов, то экспоненциальный рост станет прямой линией, двойная экспонента – простой экспонентой, степенная экспонента – обычной степенной функцией, ну а закон гиперболического роста превратится в логарифмическую зависимость. А в Excel-то как раз есть способ рассчитать все эти параметры, причем прямо на графике. Сказано – сделано. Вот что получилось:


Рисунок 3.

 

Вот они все! Причем, никакого перебора не надо, все считается по формулам математической статистики. И даже критерии Пирсона посчитаны, причем для двойной экспоненты (показана красным) совпадает с приведенным в статье. И действительно, двойная экспонента лучше всех…

Далее в статье, как и положено, следует обсуждение полученного результата. Во-первых, отмечается факт, что «гиперэкспоненциальный» рост (рост более быстрый, чем по экспоненте), лучше описывает полученные данные, чем просто экспоненциальный. Во-вторых, выражается сдержанное сожаление по поводу того, что этим «гиперэкспоненциальным» ростом не оказался рост по гиперболе, который, по словам авторов, многократно успешно применялся для аппроксимации процессов, связанных с эволюционным развитием. Таких, как рост численности человечества за последние 2 миллиона лет (!) (имеется ссылка на работы С. П. Капицы), технологического, социального и культурного прогресса, и т.п.

Почему же авторов интересует именно гиперболический рост? А потому, что в соответствующем дифференциальном уравнении, решением которого и является закон роста, скорость такого роста зависит от квадрата размера генома. Такие зависимости всегда получаются, если между элементами, число которых исследуется, имеются горизонтальные взаимные связи. Как например у нейронов в мозге: у кого нейронов в 2 раза больше, у того связей между ними в 4 раза больше. И так далее. Если где-то обнаружен гиперболический рост, то сразу ясно - это самоускоряющийся, «автокаталитический» процесс, а значит, эволюция! Но с минимальным размером генома такая штука не прошла. Закон гиперболического роста отброшен неумолимой матстатистикой.

Далее авторы оптимистично отмечают, что их наилучший закон, все-таки, не простая экспонента. При экспоненциальном законе скорость увеличения пропорциональна просто размеру генома – то есть, грубо говоря, происходит удвоение размера генома за один и тот же промежуток времени. Как, например, количество бактерий в чашке Петри. При двойном же экспоненциальном росте скорость роста пропорциональна размеру генома, умноженного на его логарифм.

Для меня экспоненциальный и гиперболический законы понятны: в первом - каждый ген со временем рождает новый, во втором - количество добавленных генов в единицу времени (например – сто миллионов лет) зависит от того, со сколькими генами каждый ген провзаимодействует. А двойной экспоненциальный? Что означает умножение размера генома на свой логарифм?

Но авторы надеются, что некая комбинация простого экспоненциального роста и продвинутого роста по гиперболе может дать полученный ими двухэкспоненциальный закон. Хотелось бы увидеть такую расчетную модель (потому что в элементарных функциях такая комбинация не выражается).

 

До сих пор излагались факты, и только факты. Дадим же теперь волю полету фантазии! Забудем о том, что наши данные отягощены неопределенностями (тем большими, чем дальше вглубь времен), систематическими погрешностями и другими факторами, которые могут вдребезги разбить любые математико-статистические построения.

У нас есть девять точек, в которых коллектив авторов во главе с крупным биологом настолько уверены, что опубликовали их в научном журнале. У нас есть санкция на их статистическую обработку с целью выявить общий закон, действующий миллиарды лет, который привел к образованию этой цепочки – ибо сами авторы этим и занимаются. И есть у нас желание найти-таки гиперболический закон роста, который бы подтвердил, что сложные системы «любят» этот закон и охотно ему следуют.

Что тут можно сделать? Посмотрим...

 

Сначала я хотел вставить в этот текст такой вот кусок:

 

Предлагаю сравнить два графика, один – повтор данных из статьи, а другой я составил сам:

 

Рисунок 4.                                                                Рисунок 5.

 

И если про законы, которыми определяется минимальный размер генома (левый график), я мало чего могу сказать, то про законы, породившие второй график, могу рассказывать долго и подробно. Построил я его, выбрав случайно (ну или почти случайно) отдельные точки из большого массива данных, полученного при имитации пуска реактора ВВЭР-1000 с помощью программы-модели. Вот как выглядит «полный» график:

Рисунок 6.

 

Глядя на «полный» график видно, что весь процесс явно разбивается на два этапа – этап ускоренного роста и этап замедленного роста (в логарифмическом масштабе). Дело в том, что, действительно, теория атомных реакторов гласит, что пока реактор в подкритическом состоянии, его мощность (при постоянной скорости ввода реактивности), растет по гиперболическому закону. После того как реактор выйдет в критику и реактивность его станет равной нулю, рост мощности будет идти по линейному закону. Красная линия, показывающая реактивность, приведена для визуального выделения этих участков.

Наличие этих двух этапов можно заподозрить, даже рассматривая выборку из редких точек на рисунке 5. При этом даже можно примерно показать пальцем, где должна быть точка перегиба.

Синяя кривая на рисунке 6 – поток нейтронов в реакторе, рассчитанный программой, учитывающей быстропротекающие процессы на мгновенных и запаздывающих нейтронах. При этом расчет приходится вести с шагом по времени порядка десятитысячных долей секунды. А есть модели попроще, но и предсказания они дают только каждая в своей области. По этим моделям выходит, что на первом участке закон будет гиперболический (к счастью, это не связано с взаимодействием нейтронов друг с другом, просто уравнение такое получается), а на втором участке – линейным. Так как параметры модели известны (подглядели в исходнике), сделаем расчет по этим моделям:

Рисунок 7

 

Вот она розовая линия – гиперболический рост по простой модели, а вот и оранжевая – линейный рост по еще более простой модели. Видно, что когда мы продляем расчет в чужую область, модели начинают врать.

Предлагаю для прикола применить к моим реакторным точкам методику, использованную в статье (прошу прощения, воспользуюсь быстрым способом, как на рисунке 3)


Рисунок 8.

 

Поразительное сходство. И коэффициенты в соответствующих формулах близкие, и критерии R2 почти совпадают, и двойная экспонента оказывается лучшим приближением. Но, к сожалению, никаких идей о настоящей природе законов, породивших эти точки, методика не дает.

Боле плодотворным представляется следующий способ: разбить точки на две группы и для каждой поискать подходящий закон. Предлагаю в первую группу включить точки с первой по пятую, а во вторую – с пятой по девятую. Так как точек в каждой группе стало меньше, то ограничимся двухпараметрическими зависимостями, которые суть – линейная, гиперболическая и экспоненциальная. Для первой группы целесообразно сразу отбросить линейный закон, а для второй – гиперболический из-за явной очевидности и для экономии времени, своего и процессорного. Попробуем...

Рисунок 9.

 

В соревновании, как и ожидалось, победили гипербола на первом участке (R2=0.9964 против 0.7144) и прямая на втором участке (1.0000 против 0.9696). Я посчитал «обобщенный» R2 этой аппроксимации (для точки под номером пять взял среднее по двум законам) – получил 0.999899. Тот же R2 для упрощенных моделей – 0.998651, ну и правильно, ведь в районе точки номер пять ни одна из них не работает.

Чем такая двухдиапазонная аппроксимация лучше, чем двойная экспонента на всем диапазоне, кроме высокого значения R2? А тем, что имея эти данные, мы можем сказать, что в исследуемом процессе в момент времени, близкий к 1469.69 секунды наступил какой-то кризис, некая критическая ситуация, в результате чего сменился закон изменения нашего параметра. Еще можем сказать, что рост параметра на втором участке определяется постоянной во времени добавкой в 1.538 миллиардных чего-то там в секунду, возможно даже, что от внешнего источника, не связанного с внутренней динамикой системы.

На самом деле: на 1445 секунде реактор достиг критики и «виртуальный оператор» прекратил ввод положительной реактивности. Дальше при нулевой реактивности (внутренняя динамика системы на нуле) рост потока нейтронов обеспечивался наличием внешнего источника нейтронов мощностью в (вы сейчас упадете!) 1.5*10-9 % номинальной нейтронной мощности реактора.  И все это по девяти точкам. Неплохо...

Поэтому, исходя из качественной близости рисунков 4 и 5, а также 3 и 8 и учитывая повышенную информативность двухдиапазонной аппроксимации простыми законами, предлагаю проделать эту же работу над данными из статьи.

(Продолжение ниже...)