Содержание к диссертации
Введение
Глава 1. Краткий обзор современных методов идентификации систем 10
1.1. Типы моделей 10
1.2. Математические модели: определение функциональной зависимости и параметризация
1.3. Приближение функций 13
1.4. Определение вида эмпирической формулы 20
1.5. Определение параметров эмпирической формулы 21
1.6. Подход к идентификации, основанный на ошибке предсказания 26
1.7. Численные методы оценивания параметров и «укороченные» методы идентификации динамических систем 31
1.8. Модели нейронных сетей 39
1.81.Классификация нейронных сетей и основные задачи, решаемые в рамках теории нейронных сетей 39
1.8.2. Модели нейронных ансамблей 41
1.8.3.Модели ассоциативной памяти и распознавания образов (сети Хопфилда) 44
1.8АМодель адаптивного порогового элемента 50
1.8.5.Нейронные сети в системах управления и идентификации динамических объектов 51
Глава 2. Использование линейных и нелинейных регрессионных моделей для задач биофизики и геобиофизики 54
2.1. Определение биологического возраста человека на основе параметров ритмической активности сердца 54
2.1.1. Параметры ритмической активности сердца как биологические маркеры старения 56
2.1.2. Линейная регрессионная модель для определения биологического возраста 64
2.1.3. Определение констант линейной модели биологического возраста и отбор значимых параметров 66
2.1.4. Квадратичная модель биологического возраста 75
2.2. Анализ палеоклиматических данных 82
2.2.1. Использование ценных палеоклиматических данных антарктической станции «Восток» для прогнозирования климата Земли. 82
2.2.2. Зависимость климатической чувствительности Земли от концентрации С02 84
Глава 3. Новый метод идентификации систем на основе критерия минимальной квадратичной невязки 88
3.1. Формулировка метода идентификации 89
3.2. Тестовые задачи идентификации для нелинейной динамической модели с предельными циклами ...94
3.3. Возможность применения метода идентификации на осноее минимальной квадратичной невязки для задач биохимии 101
3.4. Учет характера экспериментальной погрешности 105
3.5. Сравнительный анализ дифференциальных и интегральных методов
3.6. Заключение к главе 3 125
3.7. Приложение 126
Глава 4. Исследование свойств нейронных сетей как модельных динамических систем в рамках задачи идентификации 133
4.1. Основные положения 134
4.2. Нейронная сеть как статистический ансамбль 136
4.3. Модель взаимодействия популяций нейронов 141
4.4. Оценка числа стационарных состояний для нейронной сети Хопфилда 154
4.5. Задача оптимизации обучающей процедуры 158
4.6. Заключение к главе 4 160
Основные результаты и выводы 161
Список литературы 162
Список работ по теме диссертации 173
Благодарности 175
- Математические модели: определение функциональной зависимости и параметризация
- Параметры ритмической активности сердца как биологические маркеры старения
- Тестовые задачи идентификации для нелинейной динамической модели с предельными циклами
- Нейронная сеть как статистический ансамбль
Введение к работе
В настоящей работе рассматриваются различные методы идентификации и их применение для описания поведения динамических систем. Предмет теории идентификации составляет решение задачи построения математических моделей динамических систем по данным наблюдений. Идентификация системы включает ряд основных этапов: выбор математической модели на основе информации о поведении системы; получение данных наблюдений, исходя из особенностей модели; определение идентифицируемых констант модели. Теория идентификации систем является элементом общей научной методологии, что определяет актуальность данного направления исследований.
Актуальность темы. При изучении сложных нелинейных систем, таких как живой организм или отдельная клетка, практически невозможно разбить систему на простые подсистемы без потери точности ее описания. Поэтому крайне актуальной является задача разработки универсальных методов идентификации (изучения) нелинейных систем с большим числом идентифицируемых параметров.
Несмотря на широкий спектр методов описания динамических систем, существует проблема описания именно нелинейных динамических систем. Например, проведение качественных исследований систем дифференциальных уравнений позволяет определить поведение системы вблизи особой точки (положения равновесия), но не дает возможности определить параметры модели по данным наблюдений. С помощью асимптотических методов можно решать задачу идентификации лишь в линейном приближении, что часто ведет к потере точности, а в некоторых случаях, например для систем с квазистохастическими свойствами, подобная линеаризация не может быть применена в принципе без искажения основных свойств системы. Эмпирический подход применяется в тех ситуациях, когда невозможно понять и изучить внутренние Бзаимосвязн системы и поэтому он не отражает сути наблюдаемых явлений. Отметим, что существует достаточно много численных методов решения задачи идентификации нелинейных динамических систем, но их применение в сильной степени ограничено количеством идентифицируемых параметров в соответствии с вычислительными мощностями современных компьютеров.
В связи с этим перспективным представляется разработка методов идентификации нелинейных систем с большим числом параметров (порядка 10 и выше), основанных на аналитическом решении задачи идентификации. Такие методы идентификации помогут выявить устойчивые связи и закономерности, которые и будут достаточно точно описывать законы природы.
При этом не умаляется значение относительно простых линейных моделей, которые в первом приближении могут оказаться полезными при изучении ОСНОВНЫХ закономерностей, выявлении отклонений от линейных законов. Это дает возможность сделать второй шаг в изучении системы: выбрать модель из класса нелинейных моделей и определить ее параметры на основе экспериментальных данных с помощью современных методов идентификации.
Цель работы. Целью настоящей работы являлось исследование методов идентификации сложных систем на основе критерия наименьших квадратов, критерия минимальной квадратичной невязки и статистического подхода, а также их практическое применение для задач биофизики, биохимии и геобиофизики. В соответствии с целью работы были поставлены следующие задачи: построение модели биологического возраста на основе критерия наименьших квадратов, устойчивой к изменению популяционной выборки; оценка величины климатической чувствительности Земли для различных концентраций углекислого газа в атмосфере на основе аппроксимации данных наблюдений полиномами 2-й и 3-й степеней. разработка метода идентификации для широкого класса нелинейных динамических систем, позволяющего решить задачу идентификации в аналитическом виде для моделей с большим числом параметров; разработка нового статистического подхода к описанию нелинейных динамических систем, таких как нейронные сети.
Методы исследовании. Для изучения свойств математических моделей основным методом является вычислительный эксперимент. Все вычисления проводились нами с использованием встроенных функций системы MATLAB J165]. Так, для решения системы обыкновенных дифференциальных уравнений применялся метод Рунге-Кутта (функция ODE45) с использованием формул 4-го и 5-го порядков и автоматическим выбором шага на неравномерной сетке. При численном интегрировании использовалась функция TRAPZ с применением метода трапеций. Для вычисления производных вместо функции DIFF использовалась модифицированная нами схема вычисления производных в узлах разбиения, что позволило повысить точность вычислений. При вычислении обратных матриц использовалась функция (NV. Алгоритм реализации случайной шумовой компоненты основывался на использовании функции RANDN. Аппроксимация полиномами 2-й и 3-й степеней проводилась с использованием функции POLYFIT.
Наряду с численными методами исследования динамических систем использовались аналитические методы, среди которых методы идентификации систем на основе критерия наименьших квадратов, методы статистического анализа для исследования свойств нейронных сетей, в частности, Центральная предельная теорема.
Для получения параметров ритмической активности сердца использовался программно-аппаратный комплекс «Экспресс» (Центр НТТМ "Контракт") [168].
Научная новизна работы. Впервые построена линейная регрессионная модель для задачи определения биологического возраста человека, где в качестве биомаркеров старения были выбраны параметры ритмической активности сердца, которые в большой степени отражают состояние всего организма, поскольку ритмическая активность регулируется нервными и гуморальными факторами центральной нервной системы. На основе достаточно большой статистической выборки нами была получена база* данных о связи параметров ритмической активности с календарным возрастом человека, что дало возможность построить модель, устойчивую к изменению популяционной выборки. Показан принципиально нелинейный характер процессов старения, что приводит к выводу о необходимости использования нелинейных моделей с учетом того факта, что старение является сложным многопараметрическим процессом. Построена нелинейная модель биологического возраста.
На основе анализа палеоклиматических данных, полученных при бурении ледникового щита в Антарктиде, произведена оценка значения климатической чувствительности Земли для различных концентраций углекислого газа в атмосфере. Отметим, что вопрос правильного применения методов идентификации при изучении такой сложной нелинейной динамической системы, каковой является климатическая система Земли, приобретает особую значимость в связи с наблюдаемыми в последние десятилетия изменениями климата, причиной которых является парниковый эффект.
Разработан новый метод идентификации для широкого класса нелинейных динамических систем. В основе метода лежит критерий минимальной квадратичной невязки. Решение задачи идентификации получено в аналитическом виде. Метод может быть применен для моделей с большим числом параметров, что особенно актуально для исследования каскадов сложных биохимических реакций. Показана возможность применения этого метода для конкретных нелинейных динамических систем: моделей с предельными циклами (осциллятор типа Ван дер Поля и брюсселятор), модели с квазистохастическим поведением (странный аттрактор Ресслера). Показана высокая устойчивость метода к шумовой компоненте динамических параметров.
Предложено использовать видоизмененный функционал минимальной квадратичной невязки с целью учета экспериментальной погрешности измерения динамических параметров для вышеупомянутого метода идентификации, что позволяет повысить точность вычисления констант динамической модели. Эффективность такого метода демонстрируется на примере задачи идентификации констант химической реакции для системы свертывания крови.
Предложено использовать новый статистический подход, основанный на Центральной предельной теореме, для описания нелинейных динамических систем, таких как нейронные ансамбли. Подход был применен для решения ряда конкретных задач: оценки информационной емкости нейронных сетей, оценки эффективности процесса обучения и подбора параметров обучения. Статистический подход оказался также полезным для определения степени активности связанных между собой двух популяций нейронов и предсказания поведения системы в последующий момент времени.
Практическое значение работы. Разработка общих методов идентификации для нелинейных динамических систем позволяет решать специальные задачи биофизики, биологии, медицины и т. д. Например, возможно применение предложенного нами метода идентификации на основе минимальной квадратичной невязки для изучения кинетики сложных реакций (феномена колебаний рецепторного связывания или кинетики реакций системы свертывания крови), для стохастических и автоколебательных систем (реакция Белоусова-Жаботинского). Метод может найти применение для исследования многопараметрических систем, например, в теоретической медицине для разработки количественных моделей таких сложных системных заболеваний как рак, атеросклероз, артериальная гипертензия, а также для создания теории старения как совокупности патологических процессов и системных заболеваний. Возможности применения данного метода на сегодняшний день обеспечивается широким внедрением новых компьютерных и технических разработок для получения экспериментальных данных (например, ЭПР- и ЯМР- спектроскопии), позволяющих получать большие массивы информации с достаточно высокой степенью подробности.
Для нейронносетевых моделей практическое значение состоит в возможности применения полученных результатов как для понимания работы мозга, так и для создания искусственных нейронных сетей - либо на основе искусственных нейронов (нейро чипов), либо на базе компьютерных программ-имитаторов. Отмстим также, что модели нейронных сетей могут применяться для идентификации нелинейных динамических систем, путем их обучения в соответствии с набором входных и выходных сигналов.
Практическое значение разработанной нами модели биологического возраста состоит в возможности ее использования совместно с другими методами для массовых обследований населения, мониторинга окружающей среды, контроля действия лекарственных препаратов и пищевых добавок (БАДов).
Предложенный нами анализ палеоклиматических данных используется для тестирования климатических моделей. Отметим, что анализ палеоклиматических данных является практически единственным на данный момент методом проверки климатических моделей, поскольку экспериментальная проверка подобных моделей невозможна.
О содержании и структуре диссертации. В главе 1 излагаются основные методы идентификации систем. Дается анализ математических моделей, рассматривается вопрос выбора модели — определения функциональной зависимости и параметризации модели. Рассматриваются методы аппроксимации функций с использованием критерия наименьших квадратов. Обсуждается подход к идентификации, основанный на ошибке предсказания. Дается краткий анализ программных (нейронносетевых) моделей.
Несмотря на то, что набор методов идентификации систем очень широк: от простейших статистических и спектральных до методов идентификации многопараметрических нелинейных динамических систем, процедура идентификации системы является достаточно универсальной и, как правило, включает три этапа: выбор модели на основе предварительной информации о поведении системы; получение данных наблюдений, исходя из особенностей модели и экспериментальных возможностей; определение идентифицируемых констант модели и оценка степени соответствия данной модели наблюдаемым данным.
Для определения констант модели широко применяется подход минимизации функции ошибки предсказания [153], которая вычисляется в соответствии с используемым критерием и с учетом наблюдений. Подобные методы используют такие хорошо известные критерии оценки ошибки предсказания, как критерий наименьших квадратов, критерий максимального правдоподобия, критерий максимума апостериорной информации и другие.
В данной работе используется критерий наименьших квадратов для задачи определения биологического возраста и для задачи определения климатической чувствительности Земли (глава 2), а также предложенный нами критерий минимальной квадратичной невязки (глава 3).
Можно выделить в отдельную группу задачи исследования нелинейных динамических систем. В данной работе это системы, описываемые обыкновенными дифференциальными уравнениями (ОДУ) в явной форме Конщ с нелинейной правой частью (глава 3) и модели нейронных сетей, описываемые с помощью вероятностного подхода (глава 4). Асимптотические методы идентификации таких систем часто бывают неприемлемы, поскольку подобные приближения ведут к потере точности описания поведения системы. Поэтому для первой задачи (глава 3) нами предложен подход, с помощью которого константы модели находятся на уровне системы ОДУ, а для второй задачи (глава 4) предложен статистический подход на основе Центральной предельной теоремы, основным моментом которого является получение интеграла вероятностей для суммарного синаптического тока отдельного нейрона. Для таких динамических моделей путем варьирования параметров модели можно получить характерные фазовые портреты, дающие представление о поведении системы.
Отметим, что при планировании экспериментов ранее приходилось ограничиваться моделями с небольшим числом параметров. Дело в том, что среди специалистов в области идентификации систем широко распространено мнение, что «за исключением отдельных частных случаев, нелинейная задача предсказания (идентификации) не имеет конечномерных решений» [153, с. 124]. При этом имеется ввиду именно отсутствие аналитических решений. Использование же численных методов для нелинейной задачи ограничено техническими возможностями вычислительных машин.
Предложенный в работе метод идентификации на основе минимальной квадратичной невязки (глава 3) позволяет представить решение задачи идентификации широкого класса нелинейных систем в аналитическом виде. При этом задача сводится к решению системы линейных уравнений. Отметим, что для решения систем линейных уравнений существуют быстрые численные алгоритмы.
Математические модели: определение функциональной зависимости и параметризация
Выбор модели производится на основе предварительной информации о поведении системы. На практике при исследовании сложных явлений природы (физических, биологических, химических и др.) мы имеем дело с динамическими системами [166], в которых значения наблюдаемых величин зависят от их значений в предыдущие моменты времени. При этом значения параметров модели могут определяться как на уровне решения системы ОДУ (если известен вид функциональной зависимости наблюдаемых величин от времени), так и на уровне системы дифференциальных уравнений (т.е. в отсутствие информации о характере функциональной зависимости наблюдаемых величин от времени).
Вопрос о возможности определения параметров модели на уровне системы ОДУ будет рассмотрен в 3-й главе. Наиболее часто используется первый вариант, т.е. на основе полученных экспериментальным путем значений наблюдаемых величин x,(t) и информации о виде функциональной зависимости наблюдаемых величин от времени определяют конкретные константы такой зависимости- Эта информация может быть получена путем решения системы ОДУ аналитическими методами, что далеко не всегда возможно. Если аналитическое выражение функции XjQ) неизвестно, то возникает задача построения эмпирической формулы и аппроксимации.
Более подробно такая задача формулируется следующим образом. Пусть произведен ряд измерений величин х, у и в результате получен набор значений (Xf,yt) .Если аналитическое выражение функции y — f(x) неизвестно, то возникает задача: найти эмпирическую формулу у= f(x), значения которой при х = х, возможно мало отличались бы от опытных данных у,. Для большей определенности необходимо указать класс функций f(x) (линейные, степенные, показательные и т. п.). Тогда задача сводится к нахождению параметров системы. Отметим отличие данной задачи построения эмпирической формулы [175, 183] от задачи интерполирования [115, 117], которая дает значения, совпадающие в заданных точках с заданными значениями, так как график эмпирической функции у = f(x), вообще говоря, не проходит через систему точек (х Уі). Необходимо также учитывать и тот факт, что экспериментальные значения (хл у}), как правило, являются приближенными и содержат ошибки, поэтому точное совпадение эмпирической формулы в точках измерений не требуется. Кроме того, для задачи интерполирования полиномом степени т, если число узлов интерполирования и превосходит степень полинома т, то эта задача становится, вообще говоря, невозможной.
Поскольку изложенные ниже методы приближения функций, а также методы оценивания параметров эмпирических формул в большинстве случаев основаны на использовании метода наименьших квадратов (МНК), остановимся на нем подробней. Этот метод уходит корнями в классическую работу Гаусса 1809г.: «Theoria motus corporum celestium» [29]. Наиболее раннее применение МНК к временным рядам изложено в работах [91, S6].
В работе [156] доказывается, что если среди точек х0,х,,...,хпнет совпадающих и т п, то определитель системы (1.3.7) отличен от нуля, и, следовательно, эта система имеет единственное решение. Полином с такими коэффициентами будет обладать минимальным квадратичным отклонением. Метод ортогональных полиномов [9$] позволяет существенно упростить вычисления.
Коэффициенты Cj, определяемые формулой (1.3.21) представляют собой коэффициенты Фурье функции f(x) относительно заданной ортогональной системы {фДх)}[118]. Обобщенный полином с коэффициентами Фурье данной функции обладает следующими важными свойствами: 1) имеет наименьшее квадратичное отклонение от этой функции по сравнению со всеми другими обобщенными полиномами того же порядка п, 2) при увеличении числа слагаемых и младшие коэффициенты с,, как следует из формулы (1.3.21), остаются неизменными, т.е. при добавлении новых членов проделанная прежде работа сохраняется полностью, 3) при увеличении и квадратичная погрешность /„ монотонно убывает.
Подробная классификация классических алгебраических ортогональных систем (базисов) полиномов и функций непрерывного аргумента представлена в работе Дедуса Ф.Ф. и др [116]. Ортогональные базисы непрерывного аргумента представляют собой три группы базисов, объединенных порождающими их формулами и свойствами. К первой группе относятся ортогональные базисы, которые задаются общей формулой Якоби и представлены на промежутке [-1, 1]. Эта общая формула имеет параметры а и р, конкретные значения которых в различных сочетаниях приводят как к известным ортогональным системам (полиномы Лежандра, полиномы Чебышева первого и второго рода, полиномы Гегенбауэра), так и к возможным другим ортогональным полиномам. Ко второй группе относятся ортогональные базисы Сонина-Лагерра, определенные на промежутке [0, со] и имеющие всего один параметр а. Задавая различные значения параметра, можно получить множество ортогональных базисов, относящихся к группе Сонина-Лагерра. Третью группу образует ортогональный базис Эрмита, заданный на всей числовой оси (-от х со). Обычно этот базис применяется при статистических исследованиях. Эта группа представлена одним ортогональным базисом. Отметим, что классическим полиномам Якоби, Лежандра, Чебышева Сонина-Лагерра соответствуют ортогональные функции, в частности, ортогональные экспоненциальные функции,
Параметры ритмической активности сердца как биологические маркеры старения
Параметры ритмической активности сердца получают на основе данных электрокардиограммы (ЭКГ) с последующим построением графика ритмограммы, определение которой будет дано ниже.
На ЭКГ определяются [163]: изоэлекгрическая линия, - горизонтальный отрезок, записывающийся во время диастолы; зубцы — отклонения кривой вверх или вниз от изсшекірической линии (или других горизонтальных; сегментов с закруглеїшьіми или остроконечными вершинами). Различают зубец Р, отражающий г ліространение возбуждения по предсердиям, комплекс зубцов QRST, сс лветствующий возбуждению желудочков и состоящий из комплекса QRS (распространение возбуждения) и конечной части (сегмент RST и зубец Т - угасание возбуждения), а также не всегда регасгрируемый зубец U. Временные промежутки между одноименными зубцами соседних циклов носят название межцикловых интервалов или кардиоинтервалов. Ряулярность ритма определяется по межцикловым интервалам (Р-Р или R-R). Ритмограмма представляет собой упорядоченную последовательность кардиоинтервалов At, где AtRR - временной интервал между R- зубцами ЭКГ (в сек.), г номер кардиоинтервала. Ритмограмма может быть представлена в виде графика (рис. 2.1), по оси абсцисс которого откладываются номера кардиоинтервалов, а по оси ординат-временные интервалы между R- зубцами ЭКГ [124].
Процедура снятия показателей ритмической активности сердца проводилась в спокойной обстановке, в тишине, испытуемый дышал равномерно, в привычном для него темпе, не разговаривал; вставал по звуковому сигналу компьютера не резко, спокойно. Время перехода из положения лежа в положение стоя составляло 10-15 сек.
С учетом того, что модель определения БВ должна быть достаточно устойчивой к выбору состава референтной популяции, требования к состоянию здоровья испытуемых не имели жестких рамок. В обследуемую популяцию в основном входили сотрудники институтов Академии наук в подмосковном городе Пущино, - лица в возрасте от 20 до 70 лет, у которых отсутствовали выраженная сердечно-сосудистая патология и инвалидность. Объем экспериментальной выборки для мужчин и женщин составил соответственно 359 и 471 рнтмограммы.
Параметры ритмической активности были получены на программно-аппаратном комплексе «Экспресс» (Центр НТТМ "Контракт") [168]. Использовалось биполярное отведение электродов (по одному электроду на каждой руке и один электрод на правой ноге). Обработка первичных данных и вычисление параметров ритмограмм проводилась с помощью пакета программ "Экспресс" [168]. Для построения линейной регрессионной модели БВ и статистической обработки результатов использовалась выполненная автором программа на базе системы инженерных и научных расчетов MATLAB6.1.
В качестве исходного набора биомаркеров старения в нашей модели были использованы 22 параметра ритмической активности (Aj , где / — номер параметра) и один единичный параметр Ап = 1, введенный для единообразного описания свободного члена в линейной регрессионной модели (см. следующий раздел). Параметры ритмической активности сердца, используемые в модели, рассчитываются на основе ритмограммы - первичного набора кардиоинтервалов [Ai \, зарегистрированных при проведении ортостатической пробы. Часть параметров является характеристикой процессов в положении лежа (л) и стоя (с), и образует пары, например частота сердечных сокращений в положении лежа (ЧСС (л)) и стоя (ЧСС (с)). Будучи самостоятельными параметрами модели, они, тем не менее, задаются одинаковыми формулами, отличающимися только областью определения.
При спектральном анализе выявляется скрытая периодичность процессов регуляции. Согласно представлениям о многоуровневой системе регуляции сердечного ритма, развиваемой P.M. Баевским Р.М [162, 99, 100], различным диапазонам частот соответствует активность определенных уровней управления ритмом сердца. При этом выделяется два основных уровня — автономный и центральный, взаимодействующие через положительные и отрицательные обратные связи.
Параметры, An T&AIS , имеют дискретную 5-ти уровневую шкалу и представляют собой интегральные показатели, определяемые на основе довольно сложного нелинейного алгоритма распознавания с использованием фреймовых структур (пороговой логики) [168, 169]. Минимальному значению (1-й уровень) интегральных показателей соответствует нормальный резерв адаптации и функциональное состояние в пределах условной нормы, максимальному значению (5-й уровень) соответствует резкое снижение PA (Ai7 ) и ФС(Л/я ).
Данным параметрам дается следующая физиологическая интерпретация. ФС -интегральный показатель комплекса характеристик сердечно-сосудистой системы, отражающий качество работы згой системы по поддержанию гомеостаза. РЛ — интегральный показатель комплекса характеристик сердечно-сосудистой системы, характеризующий способность сердечно-сосудистой системы удовлетворять потребности организма при повышенных нагрузках.
Тестовые задачи идентификации для нелинейной динамической модели с предельными циклами
Тестирование метода производится в несколько этапов. На первом этапе получаются решения системы (3.2.1) на некотором временном интервале /є/Ч), 207 методом Рунге-Кутта с использованием формул 4-го и 5-го порядков и автоматическим выбором шага на неравномерной сетке, всего J =157 значений, которые рассматривакугся как массив данных об изменении во времени динамических параметров изучаемой системы Xj(t) (рис.3.1, а).
При вычислении обратной матрицы в формуле (3.1.10) использовалась функция INV системы MATLAB [165, стр. 145]. Восстановленные значения обозначаются а/1 2 (табл.3). Для них ошибка не превышает 2%. Имеется также возможность применения настоящего метода при наличии шума. Алгоритм реализации шумовой компоненты основан на использовании стандартной случайной функции RANDN системы MATLAB [165, стр.84]: Xl(0 = Xi(0 + yW), (3.2-4)\ где R(t) - случайная величина, распределенная по нормальному закону с математическим ожиданием 0 и среднеквадратичным отклонением 1.
Рассмотрим применение данного метода на примере тестовой задачи идентификации 3-х молекулярной системы, так называемого брюсселятора — наиболее простой модели, проявляющей свойства самоорганизации [181]. Классическим примером реальной самоорганизующейся системы, более сложной, нежели брюсселятор, является реакция Белоусова-Жаботинского. Временная зависимость динамических параметров системы уравнений брюсселятора; а -в отсутствие шумовой компоненты; б - при наличии шума 5= 0.06. Рассмотрим еще один пример использования предложенного метода для задачи идентификации системы, имеющей квазистохастичсское поведение, - аттрактора Ресслера [189]. Эта задача интересна тем, что необходимо идентифицировать поведение системы, для которой отсутствует предельный цикл или устойчивый фокус в фазовом пространстве, но траектория движения точки занимает ограниченный объем этого пространства.
Для задач идентификации систем с предельными циклами в целях упрощения модели часто используется метод линейного приближения для динамических уравнений, что значительно упрощает подбор параметров модели. Например, при исследовании динамики реакции Белоусова-Жаботинского рассматриваются только линейные члены динамического уравнения [181, с. 127-130 ].
Для квазистохастических систем метод линейного приближения динамических уравнений неприменим, вследствие принципиального отсутствия у подобных систем линейного асимптотического поведения. В отличие от других, наш метод позволяет использовать общий вид нелинейной динамической модели (3.1.1) и определять константы динамической модели в аналитическом виде.
При изучении сложных нелинейных динамических систем обычно используются подход, основанный на разделении сложной системы на простые подсистемы. Подобная схема работает лишь в тех случаях, когда возможно экспериментальное выделение фрагментов исходной системы без необратимого изменения их свойств. Однако такая возможность имеется далеко не всегда.
Например, при изучении зависимости скорости ферментативных реакций от концентрации субстрата часто наблюдаются отклонения от гиперболического закона [36] , причиной которых может быть достаточно сложный механизм ферментативного процесса: неэквивалентность субстратсвязывающих центров или взаимодействия между субстратсвязывающими центрами в молекуле олигомерного фермента, наличие в ферментном препарате форм фермента (в том числе олигомерных форм), различающихся по каталитическим свойствам [57]. К тому же многие биохимические реакции сопровождаются образованием целого спектра промежуточных продуктов.
Довольно часто экспериментатору бывает трудно сделать окончательный выбор среди возможных механизмов, предлагаемых для объяснения наблюдаемых кинетических аномалий. Поэтому остаются популярными эмпирические уравнения, предлагаемые для описания кинетики ферментативных реакций [1,152]. Как мы видим, даже при изучении отдельных ферментативных реакций исследователи сталкиваются с нетривиальной задачей выделения простых подсистем. Когда же речь идет о каскадах биохимических реакций, сопровождающих, например, такие явления, как аноптоз [174] или свертывание крови [45] , задача выделения простых подсистем становится еще более нетривиальной. Предлагаемый новый метод идентификации систем позволяет строить математические модели сложных (многопараметрических) нелинейных явлений (а именно с такими явлениями, как правило, приходится иметь дело в биологических системах) без экспериментального исследования отдельных фрагментов исходной системы.
В биохимии для анализа кинетики ферментативных реакций широко применяется ряд методов идентификации, среди которых следует отметить графический метод [108, с. 695-699] и метод Проии [155], как наиболее распространенные. Не углубляясь в детальное описание, отметим, что эти методы, а также и другие [108, с. 372-375], которые используются в настоящее время для изучения кинетики ферментативных реакций, позволяют работать только с линейными или асимптотически линейными системами. Однако подавляющее большинство ферментативных реакций описывается системами нелинейных дифференциальных уравнений. Это означает, что приходится либо специально подбирать условия эксперимента, чтобы иметь возможность линеаризовать систему, например, значительно увеличивая концентрацию одного из реагирующих веществ, либо изучать реакции с сильно отличающимися скоростями протекания. В последнем случае, в силу теоремы Тихонова [108, с. 707-708], удается уменьшить число степеней свободы идентифицируемой системы.
Как уже отмечалось, для биохимических задач константы описанной выше динамической модели aflJl можно интерпретировать как константы биохимических реакций, а динамические параметры изучаемого явления Xft)- как концентрации участвующих в реакциях веществ, включая концентрации промежуточных продуктов реакции.
Рассмотрим более подробно применение предложенного метода к описанию конкретной биохимической реакции лиганд-рецепторного взаимодействия. Поскольку описание модели простого рецепторного связывания не вызывает трудностей [108, с. 342], представляется интересным рассмотрение более сложной модели лиганд-рецепторного взаимодействия, учитывающей процессы деградации лигандов (разрушения ферментами) и интернализации (экзоцитоза) рецепторов, лигандов или лиганд-рецепторных комплексов.
Нейронная сеть как статистический ансамбль
Рассмотрим нейронный ансамбль - совокупность взаимодействующих пороговых элементов, динамика которых описьшается уравнением эволюции (4.1.1). На вход каждого нейрона подается внешний квазистохастический бинарный сигнал, который принимает значения «О» или «1» в дискретные моменты времени.
Представление нейронной сети как статистического ансамбля основывается на Центральной предельной теореме. В классическом варианте Центральной предельной теоремы суммы от большого числа независимых или слабо зависимых случайных величин имеют распределение вероятностей, близкие к нормальному распределению, если выполняется условие конечности математических ожиданий и дисперсий независимых случайных величин.
Таким образом, мы получили еще одно важное соотношение для дисперсии суммарного синаптического тока, которое вместе с выражением для математического ожидания б (4.2.2), будет использовано для дальнейших рассуждений.
Вероятность того, что нейрон может перейти из состояния покоя (St-O) в активное состояние (5Ї=1), Р , выражается интегралом от произведения функции распределения (4.2.8) для Д9 и функции распределения (4.2.1) для G . Результаты математического моделирования для эволюции популяции нейронов, описываемой уравнением (4.3.2), для различных значений параметра Dk при фиксированных значениях остальных параметров: 9, =100, N = 20000, представлены на рис. 4.3 а, 6, в, г. Параметр Мк подбирался таким образом, чтобы обеспечить равновесное состояние в точке A(t) = \f2. Для данных значений BnOjt,=100 и N = 20000 и Мк =0.01 функция интеграла вероятностей erfc(x) = l в этом случае активность популяции нейронов в последующий момент времени равна A(t + є) = 112. На рис. 4.3, а видно, что при различных начальных условиях 0 A(tu)S 1, при Dk =3.5, имеется одно стационарное состояние А = 0.5. При уменьшении значения параметра Z)t=1.2 появляется еще одно стационарное состояние А 0. При Dk =0.5 имеется одно состояние неустойчивого равновесия А =0.5 и два стационарных состояния А = 0.62 и А = 0. При Dk=0.l остается состояние неустойчивого равновесия А = 0.5 и стационарное состояние при А = 0, появляется стационарное состояние Л = 1.
Посредством варьирования параметров уравнения эволюции (4.3.8,4.3.7) можно получить характерные картины динамического поведения нейронных ансамблей. Сначала подбираются параметры для случая, когда воздействие внешнего источника отсутствует MS]a=07 DS!a=Q.
В отсутствие внешнего сигнала существует стабильный предельный цикл (рис. 4.4, а). При наличии достаточно сильного внешнего сигнала колебательная активность исчезает (рис. 4.4, б, в). Наиболее ярким примером такой активности нейронов в мозге человека может быть а-ритм. По данным ЭЭГ различают несколько типов электрической активности мозга, из них в состоянии бодрствования наиболее ярко выражен а-ритм (для здорового человека его частота -10 Гц). Он регистрируется на ЭЭГ в затылочной области, соответствующей зоне зрительной коры, когда у испытуемого закрыты глаза. При открытых глазах, когда на зону зрительной коры по зрительным волокнам поступает внешний сигнал, а-ритм исчезает. Кроме того, для а-ритма весьма характерной является картина, когда этот ритм представляет собой амшштудно модулированный сигнал, так называемые «веретена» а-ритма. Если испытуемый смотрит на периодические вспышки света, то регистрируется а-ритм более высокой амплитуды в том случае, когда частота вспышек совпадает с частотой а-ритма испытуемого. Таким образом достигается явление резонанса для этого типа ритмической активности мозга. Ниже представлены результаты моделирования такого сигнала.
Если допустить, что внешний сигнал меняется во времени, например по гармоническому закону, то для случая с предельным циклом можно получить картину изменения во времени амплитуды предельного цикла. На рис. 4.5, а показан фазовый портрет для системы из двух взаимодействующих популяций нейронов для случая с предельным циклом, для которого статистические параметры весовых коэффициентов, количество нейронов и значение Q„op остаются теми же, что и для предыдущего случая с предельным циклом (рис.4.4, а), а внешний сигнал является переменным MS}2-A-cos(2-n-v), где А = 4.5, v = 0.0055, интервал времени эволюции популяций t=[0, 600]. Частота внешнего сигнала подобрана таким образом, чтобы в один период изменения амплитуды внешнего сигнала укладывалось несколько периодов колебательной активности ансамбля из двух взаимодействующих популяций (в данном случае 10). В начальный момент времени активности первой и второй популяции имеют значения А1 = 0, Аг = 0.5. В последующие моменты времени наблюдается колебательная активность каждой из популяций. Из рисунка видно, что для такого достаточно низкого значения параметра амплитуды А внешнего сигнала предельный цикл сохраняется, не переходя в фокус, но амплитуда предельного цикла меняется периодическим образом. Если рассмотреть картину изменения во времени активности одной из популяций, например первой (рис 4.5, б), то можно видеть амплитудно-модулированный сигнал. Такая картина амплитудно-модулированного сигнала характерна для электрической активности мозга, например, для диапазона частот в области а - ритма.
