Прикладные задачи
нелинейной теории колеоании и волн
УДК 530.182 https://doi.org/10.18500/0869-6632-2020-28-4-397-413
Реконструкция уравнений нейроподобного осциллятора, моделируемого системой фазовой автоподстройки частоты с запаздыванием, по скалярному временному ряду
М.В. Сысоева1,2, И. В. Сысоев1&3&4, В. И. Пономаренко1&3, М.Д. Прохоров1
Россия, 410012 Саратов, ул. Астраханская, 83 4Национальный исследовательский Нижегородский государственный университет имени Н.И. Лобачевского
Россия, 603950 Нижний Новгород, пр. Гагарина, 23 E-mail: bobrichek@mail.ru, ivssci@gmail.com, ponomarenkovi@gmail.com, mdprokhorov@yandex.ru Автор для переписки Сысоева Марина Вячеславовна, bobrichek@mail.ru Поступила в редакцию 22.04.2020, принята к публикации 8.06.2020, опубликована 31.08.2020
Цель настоящего исследования - разработка методики реконструкции уравнений нейроподобного осциллятора, описываемого моделью системы фазовой автоподстройки частоты с запаздыванием, по скалярному временному ряду наблюдаемой. Методы. Располагая скалярным рядом только одной переменной, соответствующей трансмембранному потенциалу, для восстановления вектора состояния ещё одна переменная получается численным дифференцированием со сглаживанием полиномом, а третья - численным интегрированием методом Симпсона. Далее вводится целевая функция, описывающая длину нелинейной функции при пробном времени запаздывания, и проводится её минимизация. Результаты. Предложенным методом удаётся реконструировать время запаздывания, эффективные параметры системы и нелинейную функцию модели в различных периодических и хаотических режимах, включая режимы перемежаемости. Метод работоспособен в том числе при наличии 1-процентного измерительного шума. Заключение. Описанный метод полезен как средство реконструкции моделей нейронов по экспериментальным данным внеклеточных или внутриклеточных записей из мозга или в культуре.
Образец цитирования: Сысоева М.В., Сысоев И.В., Пономаренко В.И., Прохоров М.Д. Реконструкция уравнений нейроподобного осциллятора, моделируемого системой фазовой автоподстройки частоты с запаздыванием, по скалярному временному ряду//Известия вузов. ПНД. 2020. T. 28, № 4. С. 397-413. https://doi.org/10.18500/0869-6632-2020-28-4-397-413
Статья опубликована на условиях лицензии Creative Commons Attribution License (CC-BY 4.0).
Финансовая поддержка. Работа выполнена при поддержке Российского фонда фундаментальных исследований, грант № 19-02-00071 (реконструкция уравнений нейроподобного осциллятора) и в рамках государственного задания (разработка новых систем с хаотической динамикой).
© Сысоева М.В., Сысоев И.В., Пономаренко В.И., Прохоров М.Д., 2020
https://doi.org/10.18500/0869-6632-2020-28-4-397-413
Reconstructing the neuron-like oscillator equations modeled by a phase-locked system with delay from scalar time series
M. V. Sysoeva1&2, I. V. Sysoev1&3&4, V.I. Ponomarenko1&3, M.D. Prokhorov1
The purpose of this work is to develop the reconstruction technique for the neuron-like oscillator equations descibed by a phase-locked system model with delay from scalar time series. Methods. We reconstruct the state vector given a scalar series of only one variable corresponding to the transmembrane potential. The second variable is obtained by numerical differentiation with smoothing by a polynomial. The third variable is obtained by numerical integration using the Simpson method. Next, the target function describing the length of the nonlinear function at a trial delay time is constructed and minimized. Results. The delay time, effective system parameters, and nonlinear function can be reconstructed using the proposed method. The method gives correct results in various periodic and chaotic regimes, including intermittency. The method works even in presence of 1% measurement noise. Conclusion. The described method is useful as a tools for reconstructing neuron models from experimental data of extracellular or intracellular recordings from the brain or in culture.
Reference: Sysoeva M.V., Sysoev I.V., Ponomarenko V.I., Prokhorov M.D. Reconstructing the neuron-like oscillator equations modeled by a phase-locked system with delay from scalar time series. Izvestiya VUZ. Applied Nonlinear Dynamics, 2020, vol. 28, no. 4, pp. 397-413. https://doi.org/10.18500/0869-6632-2020-28-4-397-413
This is an open access article distributed under the terms of Creative Commons Attribution License (CC-BY 4.0).
Acknowledgements. This work was supported by the Russian Foundation for Basic Research, grant No. 19-02-00071 (reconstruction of equations of neuron-like oscillator) and within the framework of the state task (construction of new systems with chaotic dynamics).
Введение
Идея реконструкции математических моделей по временным рядам динамики описываемых ими объектов достаточно давно стала весьма привлекательной [1]. Если реконструкция удастся, можно решить большое число важных задач: верифицировать существующие модели, измерить недоступные прямому измерению величины (как переменные, так и параметры), восстановить архитектуру связей в сложных сетях колебательных элементов, сделать прогноз будущего поведения, произвести кластеризацию объектов и многое другое [2]. Основная проблема здесь в том, что задача реконструкции - классическая обратная некорректная задача, для решения которой, как правило, не хватает априорных знаний и достаточного объёма измерений.
Попытки решить задачу о «чёрном ящике» - восстановить уравнения в самом общем виде только на основе измеренных сигналов, уже давно были признаны бесперспективными [3] и все основные результаты были получены для некоторых узких классов систем [4,5]. В то же время, излишне высокие требования к априорной информации, часто характерные для методов работы
со скрытыми переменными [6-8], также могут навредить при практическом приложении методов к экспериментальным данным, поскольку реальные системы никогда не описываются моделью в точности, а метод может стать излишне «ломким» - даже небольшие отклонения оператора эволюции от предполагаемого приведут к невозможности реконструкции по временным рядам.
Ещё одна типичная проблема при реконструкции - недоступность части переменных модели: как правило, для реконструкции доступна единственная наблюдаемая. Если оператор эволюции предполагается в самом общем виде (просто набор дифференциальных или разностных уравнений), эту проблему принято решать методом временных задержек или методом численного дифференцирования [9]. Но если имеется задача построить адекватную модель с осмысленными параметрами, приходится прибегать к более сложным в техническом отношении подходам, в частности, методам реконструкции при наличии скрытых переменных [6]. Богатство динамики при возможности получить все необходимые временные ряды по одной наблюдаемой также привело к большому интересу к системам с запаздыванием, для которых за последнее время было разработано большое число специализированных подходов [10-18], в том числе при наличии запаздывания в связях [19,20].
Если восстанавливаемые уравнения динамики содержат неизвестные нелинейные функции, в том числе функции связи [21], их часто раскладывают в функциональные ряды, например, в полиномиальные [22] или тригонометрические [23]. Такой подход имеет ограничения, поскольку увеличение числа свободных неизвестных параметров и не самые хорошие свойства используемых функций, например, их неограниченный рост на краях измерительного диапазона, ведут к ошибкам и повышению требований к длине рядов и частоте выборки. Частично решить эту проблему можно, перейдя к методам, основанным на табличной реконструкции хотя бы одной нелинейной функции за счёт специального выбора целевой функции, используемой в качестве критерия точности подгонки параметров, как это было предложено в [24].
Таким образом, для практического применения очень важным является выбор модели. В частности, для реконструкции моделей нейронов по экспериментальным данным на наш взгляд весьма перспективна система, предложенная в работах [26-28], поскольку она воспроизводит некоторые основные типы нейронной активности, по одной наблюдаемой можно реконструировать все три переменные, при этом одну из них - методом численного дифференцирования, а другую - интегрирования, что должно понизить чувствительность к шумам измерений по сравнению с вариантом, когда обе переменные пришлось бы реконструировать дифференцированием. Если построить целевую функцию как в работах [21,24] на основе длины описания нелинейной функции, проблема с тем, что при интегрировании результат будет определён с точностью до константы, снимается, а метод приобретёт дополнительную общность на случай, если в реальной системе нелинейная функция не совсем такая, как в модели. Количество описываемых моделью режимов можно существенно повысить, если ввести в неё небольшое запаздывание, как показано в [29]. При этом введение запаздывания не создаёт дополнительных проблем с реконструкцией вектора состояния.
Целью данной работы является разработка и апробация в численном эксперименте методики реконструкции по скалярному временному ряду времени запаздывания, нелинейной функции, параметров и скрытых переменных предложенной в [29] модели.
колебательной активности, качественно отображающие некоторые режимы изменения мембранного потенциала нейрона, например, регулярную импульсную активность и пачечные разряды с различным числом импульсов в пачке, а также квазихаотические колебания [26]. В работе [27] было произведено разбиение плоскостей параметров модели на области существования динамических режимов с различным количеством импульсов в пачке, были исследованы бифуркационные механизмы переходов между режимами модели. В настоящей работе предлагается использовать модель системы фазовой автоподстройки, модифицированную путём введения в неё зависимости от запаздывающего аргумента. Исследованию влияния запаздывания на динамику системы ФАП посвящено множество работ, в частности [30,31]. Введение запаздывания позволяет дополнительно существенно обогатить динамику системы, в которой появляются различные примеры хаотического перемежающегося поведения, как это показано в [29]. Имеется довольно большое число работ, в которых при моделировании ионных токов и динамики электрической активности нейронов используются уравнения с запаздыванием [32-35]. В этих работах запаздывание вводится в модели, изначально записанные из биофизических принципов. Рассматриваемая нами модель, предложенная в [29], является феноменологической. Она имеет некоторое сходство с моделями, записанными из первых принципов, по структуре уравнений, но получена из радиотехнической модели ФАП [27] путём добавления запаздывания в одну из переменных, которая может рассматриваться как трансмембранный потенциал [28]. Именно эта переменная по аналогии с биофизическими моделями является измеряемой и она же может быть использована для связи между генераторами (нейронами), поэтому в [29] запаздывание было добавлено именно в неё.
Динамика рассматриваемой модели описывается следующей системой дифференциальных уравнений:
Tt = »•
% = - О
£1£2 ^ = Y - (¿1 + £2)z - (1 + £1 cos ф)y(t - т),
где ф - текущая разность фаз подстраиваемого и опорного генератора, у - начальная частотная расстройка, £1,^2 - параметры инерционности фильтров. Применительно к динамике нейрона переменную у можно интерпретировать как описывающую изменение мембранного потенциала, параметры £1 и £2 позволяют задавать необходимый динамический режим, у оказывает воздействие, сходное с воздействием внешнего тока в модели Ходжкина-Хаксли [36], а т соответствует времени рефрактерности нейрона после генерации потенциала действия.
На рис. 1 построена бифуркационная диаграмма зависимости переменной у от времени запаздывания т. График построен сечением фазового пространства системы ( ) плоскостью z = 0 при параметрах, соответствующих режиму пачечных разрядов с малым числом импульсов в пачке (у = 0.075, £1 = 4.5, £2 = 10). Видно, что с изменением параметра т модель демонстрирует как регулярные колебания разного периода, так и нерегулярные колебания. С увеличением времени запаздывания в системе появляется мультистабильность (или, по крайней мере бистабильность). Так, одновременно со сложными хаотическими режимами, включая режимы перемежаемости, при значениях 4.125<т<4.875 существует режим простых, близких к гармоническим колебаний.
Основываясь на этой бифуркационной диаграмме, для дальнейшего тестирования предложенного метода реконструкции использовали периодический режим при т = 2, слабохаотический режим при т = 2.71875, режим развитого хаоса при т = 3.125 и режим с перемежаемостью
Рис. 1. Бифуркационная диаграмма системы ( ) по параметру т при фиксированных значениях параметров у=0.075, ti = 4.5, £2 = 10, которые соответствуют режиму сильно нелинейных периодических колебаний при т = 0. Диаграмма построена с использованием сечения пространства состояний плоскостью г = 0, отложена переменная у
Fig. 1. Bifurcation diagram for parameter x in the system ( ). Values of parameters y = 0.075, £i = 4.5, £2 = 10 are fixed and correspond to the regime of strongly nonlinear periodic oscillations at x = 0. The diagram is constructed using a cross section of the state space with the plane г = 0, the variable y is shown
при т = 4.78125. При больших значениях т система демонстрирует динамику, которая определяется почти исключительно величиной т и не наследует свойств оригинальной системы без запаздывания [27].
£1 +g2 tlt2
где /(ф) - неизвестная в общем случае нелинейная непрерывная функция. Таким образом, то, что ф получается численным интегрированием с точностью до константы, не имеет значения. При этом далее будем рассматривать свёрнутую фазу ф такую, что 0 ^ ф < 2п.
Следуя идеям работы [24], будем строить целевую функцию для вычисления коэффициентов ао и а1, исходя из минимизации длины описания нелинейной функции /. Это позволит отказаться от явного разложения функции / в ряд и, таким образом, одновременно повысит общность метода (он может использоваться для произвольной /) и уменьшит параметризацию (число подлежащих оценке параметров а), что в свою очередь улучшит статистические свойства оценок оставшихся в модели коэффициентов (чем меньше параметров, тем лучше их оценки). Перепишем уравнение (2), выразив из него /:
О / \\ 1 ^ 1 ^ ^ \\
/(ф) =аоЙТ-Г) +а1 - ЙГ^та- (3)
Временной ряд (Ы/(И также рассчитаем численным дифференцированием.
Далее введём сортирующее отображение ^(п), где п - номер точки в ряде, ставящее в соответствие п-й точке в исходном ряде её номер Q(п) в ряде, отсортированном по возрастанию ф. Рассмотрим также обратное отображение Q-1, вычисляющее по номеру в отсортированном ряде номер в исходном, так что Q-1(Q(n)) = п. Рассмотрим точку, находящуюся в отсортированном ряде непосредственно перед Q(n)-й, тогда в исходном ряде она имеет номер Q-1(Q(n) — 1), который обозначим как рп для краткости. Тогда приращение функции / на отрезке [ф(рга); ф(п)] будет выражаться формулою ( ), где 6 = т/ (Д£):
Д^(п) = 6) I(п) — шЬё) §ы, (4)
Ду-1 (п — 6)= 1 1
Ди(п) =
у(п — 6) у(рп — 6) & г(п) г(рп)
у(п — 6) у(рп — 6)&
При этом значение п такое, что Q(n) = 0 (будем здесь и далее считать, что нумерация начинается с нуля, то есть п = 0,1,..., N — 1), не позволено, поскольку для него нет соответствующего рп (нет предыдущей точки в отсортированном ряде — наше значение самое маленькое). Тогда в качестве целевой функции рассмотрим величину Ь
L(ао, а0 = ^ 8п- (5)
Из (4) очевидно, что Ь будет зависеть от ао и а1 квадратично. Значит, формулу (5) можно рассматривать как формулировку задачи на наименьшие квадраты для аппроксимации величин ДЦп) величинами Ду-1(п — 6) и Ди(п), которые можно рассматривать как базисные функции. При этом Ьп станут соответствовать невязкам. При верном выборе ао и а1 целевая функция Ь будет много меньше, чем при неверном, то есть когда параметры не угаданы и функция / имеет разрывы почти в каждой точке.
Поскольку даже при верном выборе ао и а1 в общем случае L > 0, то при конечном N полученные таким образом оценки должны быть смещёнными - их математическое ожидание не будет в точности совпадать с истинными значениями. Если функция / непрерывна, то Sra^0 при N ^ те, поскольку 0 ^ ф < 2п и, значит, бесконечное число измеренных значений ф(п) будут плотно заполнять конечный полуинтервал. По сути, если дополнительно считать, что функция не только непрерывна, но и дифференцируема на полуинтервале 0 ^ ф < 2п, то Sra ^ df (ф(рп)), то есть Ьп является аппроксимацией дифференциала функции / в точке ф(рп) справа (или в точке ф(п) слева, что для дифференцируемой функции в пределе одно и то же). Тогда lim^L =
р 2 р 2 2 2 = J (df) = J (4/У^ф) ( ^ф) = (df/dф) dф ^ 0, то есть оценки предложенным методом являются асимптотически несмещёнными.
Предложенный подход позволяет оценить значения ао и а1 , являющиеся комбинациями исходных параметров у, 81 и 82, но не сами параметры. Следует отметить, что, если рассмотреть третье уравнение системы ( ), то можно видеть, что при замене 1 + 81 cos ф = /(ф) оценить отдельно у, 81 и 82 будет нельзя, потому что все члены уравнения, включая правую часть, стоят с неизвестными коэффициентами (или функциями). Таким образом, система уравнений для определения коэффициентов, которую можно было бы составить, оказалась бы вырождена (нет свободного члена). Поэтому представленный метод ничем не уступает подходу [22], где функция (ф) аппроксимировалась бы напрямую.
В приложении к системе (1) отрезок, на котором перебирается пробное запаздывание т, должен начинаться от 0 и продолжаться до значений, соответствующих таким т, при которых динамика заведомо не похожа на наблюдаемую. Если рассматривать режимы, наследующие основные свойства системы без запаздывания, имеет смысл рассматривать достаточно малые т, как это показано в [29], в частности, т < 6 для режимов, рассмотренных в работах [27,28].
I llillffllii Iii
О 500 1000 1500 2000 2500
Рис. 2. Временной ряд наблюдаемой у и восстановленные временные ряды переменных q> и z при т = 3.125 Fig. 2. The time series of the observed variable у and the reconstructed time series of the variables ф and г for x = 3.125
Были выбраны 4 режима при различных т со следующими характеристиками:
Во всех рассмотренных случаях истинные значения параметров были одни и те же: у = 0.075, t~i = 4.5,12 = 10, что соответствует а0 = 1/600 = 0.001(6) иа1 = -29/90 = -0.3(2).
Визуально временные ряды переменных z и ф (рис. ), реконструированные методами численного интегрирования и дифференцирования соответственно, похожи на исходные.
Восстановленные при найденном т параметры (ао и ai) и значение целевой функции ( ) представлены следующим списком:
На рис. , b изображена истинная нелинейная функция (серая линия) и восстановленная нелинейная функция (чёрные точки). Неидеальное совпадение графиков истинной и восстановленной нелинейных функций обусловлено тем, что численное интегрирование производится с точностью до константы - графики могут быть смещены по горизонтали на произвольное число.
Восстановленная нелинейная функция (рис. , Ь) имеет «выбросы» - выбивающиеся, сильно смещённые оценки нелинейной функции, которые получаются при y(t — т) близких к нулю.
Рис. 3. Восстановление параметров модели без шума, а - Зависимость значения целевой функции от пробного времени запаздывания £(т). Минимум зависимости соответствует восстановленному времени запаздывания, равному истинному т = т = 3.125 и обозначено чёрной звёздочкой; Ъ - Восстановленная нелинейная функция (серые точки) и истинная нелинейная функция (чёрная линия)
Fig. 3. Reconstructed model parameters in the absence of noise, a - Dependence of the target function on the trial delay time L(x). The minimum of the target function is observed at the reconstructed delay time equal to the original one x = x = 3.125 and is indicated by a black asterisk, b - The reconstructed nonlinear function (gray dots) and the original nonlinear function (black line)
Это вызвано тем, что в формуле ( ) при у(п — 6) = 0 или при у(рп — 6) = 0 появляется деление на ноль. Хотя строго значение 0 в наблюдаемых рядах отсутствует, любые очень малые у фактически опасны для процедуры, так как матрица базисных функций, используемая при решении задачи на наименьшие квадраты, становится плохо определенной. Особенно это должно проявляться при наличии шума, так как не равные нулю значения у могут стать близкими к нулю из-за шума. Поэтому ещё одним важным параметром является степень ограничения на малость значений переменных в знаменателях в формуле ( ). Назовём этот параметр ц, он будет зависеть от дисперсии сигнала. Все слагаемые Ъп с такими п. для которых \\у(п — 0) < и или Iу(Рп — б)| < М- исключим из суммы в формуле ( ), используемой для целевой функции L. Малое (х, очевидно, неприемлемо: сохранится слишком много «плохих» значений. Если взять очень большое (х, то можно уменьшить ошибки, связанные с возможным делением на малое число, но при этом пожертвовать частью данных и, таким образом, ухудшить статистические свойства оценок параметров и нелинейной функции. В итоге, экспериментально было определено, что значение |х = 0.2 является оптимальным, поскольку при этом выбрасывается менее половины значений исходного ряда (длина ряда уменьшается некритично), но все у(п — 0) оказываются примерно одинаковыми по порядку величины.
• о \\ 1
Рис. 4. Зависимость ошибки аппроксимации от количества точек т, используемых для сглаживания при дифференцировании, при различных режимах, отличающихся параметром т: т = 2 (а), т = 2.71875 (Ь), т = 3.125 (с), х = 4.78125 (d)
Fig. 4. The dependence of the prediction error on the number of points used for smoothing at differentiation in various regimes that differ by the parameter т: x = 2 (a), x = 2.71875 (b), x = 3.125 (с), x = 4.78125 (d)
что может быть объяснено возникновением новых, более быстрых временных масштабов в сигнале с увеличением времени запаздывания.
Полученные таким образом значения пг позволяют верно реконструировать время запаздывания для всех рассмотренных режимов. Необходимое количество пг и восстановленное при этом значении время запаздывания т для каждого режима представлены следующим списком:
На рис. 5 приведены результаты оценки времени запаздывания (рис. : , а) и реконструированная нелинейная функция (рис. 5, Ь) при наличии в ряде наблюдаемой 1-процентного измерительного шума. Показаны результаты при т = 3.125. Для построения кривой Ь(т) использовалось значение пг = 125, подобранное путём построения зависимости Ь(пг), как это показано на рис. , с. Видно, что зависимость Ь(т) сохраняет при наличии шума тот же качественный вид, что и в его отсутствии (см. рис. , а), но минимум становится менее выражен, значение целевой функции в минимуме больше, а положение минимума оказывается смещено
Рис. 5. Результаты реконструкции с 1-процентным измерительным шумом, а - Зависимость значения целевой функции от пробного времени запаздывания £(х). Минимум зависимости соответствует восстановленному времени запаздывания т = 3.15625, истинное время запаздывания в этом случае х = 3.125 и обозначено чёрной звёздочкой. Ъ - Восстановленная нелинейная функция (серые точки) и истинная нелинейная функция (чёрная линия)
Fig. 5. Reconstruction results in the presence of 1% measurement noise, a - Dependence of the target function on the trial delay time L(x). The minimum of the target function is observed at the reconstructed delay time equal to the original one x = x = 3.125 and is indicated by a black asterisk, b - The reconstructed nonlinear function (gray dots) and the original nonlinear function (black line)
вправо на один шаг интегрирования At. Аналогичные различия между зависимостями целевой функции от пробного времени запаздывания были ранее отмечены в ансамбле систем первого порядка в работе [38]. Реконструированная при наличии шума нелинейная функция представляет собою как бы несколько наложенных кривых, различным образом примерно воспроизводящих оригинальную нелинейную функцию - см. рис. f, Ъ. Такое поведение, по-видимому, является следствием хаотической природы колебаний, когда одним и тем же значениям переменной ср могут соответствовать различные по форме колебательные паттерны, и, как следствие, различные значения у и г. Поскольку такие паттерны имеют разный характерный временной масштаб и форму, сглаживание при численном дифференцировании будет вносить для каждого из них разные искажения.
Заключение
В данной работе на основе уже сформулированных ранее идей нами предложен метод реконструкции модели системы фазовой автоподстройки с запаздыванием в собственной динамике [29], основанной на работах [26-28,37]. Данная система интересна в первую очередь уникальным сочетанием богатства динамических режимов и удобства для реконструкции. Ранее было показано, что даже для систем ФАП с фильтром первого порядка введение времени запаздывания приводит к возникновению дополнительного временного масштаба колебаний [39,40]. В случае модели, предложенной в [29], введение запаздывания дало возможность помимо режимов генерации отдельных импульсов (спайков) и пачек импульсов (беретов) добиться также наличия режимов перемежаемости.
Предложенный метод использует табличную реконструкцию нелинейной функции системы, причём в качестве целевой функции для подгонки параметров модели выступает величина, связанная с длиною нелинейной функции. Временные ряды недостающих переменных получаются численным интегрированием и дифференцированием. При этом процедура реконструкции не использует решение уравнений динамики и не предполагает стартовых догадок о значениях параметров, как это делается в методах, основанных на работе со скрытыми переменными [6,8].
Поэтому её успех зависит исключительно от свойств наблюдаемой траектории: её регулярности и зашумлённости. Нерегулярные режимы реконструируются лучше, а регулярные - хуже, но их сосуществование (мультистабильность) не влияет на качество реконструкции - аналогичная ситуация имеет место и для режимов, имеющих существенно различные времена запаздывания.
С использованием предложенного подхода удаётся реконструировать время запаздывания, эффективные параметры и нелинейную функцию модели, в том числе при наличии измерительного шума. За счёт табличного подхода к реконструкции нелинейной функции удаётся нивелировать ошибку на константу при определении одной из переменных модели путём численного интегрирования (соответствующая константа сокращается при построении целевой функции), а также расширить класс рассматриваемых систем: в действительности, если нелинейная функция системы отличается от рассматриваемой, метод всё равно будет работать, так как требует только предположений о её непрерывности.
В отличие от ранее разработанных подходов для систем с запаздыванием [24,38] и ней-роосцилляторов с запаздывающими связями [19,20], данный метод позволяет реконструировать систему третьего порядка по одной скалярной переменной. При этом чувствительность метода к шумам не хуже, чем подхода, ранее разработанного для систем второго порядка без запаздывания - обобщённых осцилляторов ван дер Поля [41]. Это достигается за счёт того, что только одна из переменных модели реконструируется по наблюдаемой методом численного дифференцирования, а другая - методом численного интегрирования. При этом впервые в подходе, основанном на построении целевой функции задачи оптимизации на основе длины описания нелинейной функции системы, рассматривается нелинейная функция от ненаблюдаемой переменной.
Поскольку при реконструкции временные ряды {zn} и {dz/dt}n приходится получать методом численного дифференцирования, сглаживание исходного ряда оказывается критически важным. Ранее выбор количества точек для сглаживания также обсуждался, например, в работах [20, 38], но чёткие критерии не были сформулированы. В данной работе нам удалось показать, что построенная нами целевая функция, изначально предназначенная для расчёта коэффициентов модели и табличной реконструкции нелинейной функции, годится не только для определения времени запаздывания, как это уже было в [38], но и для определения оптимального числа точек, по которым производится сглаживание. Такой подход потребовал значительного увеличения времени счёта, но при использовании современной вычислительной техники и возможности распараллелить расчёты, это не составляет большой проблемы.
Предложенный нами метод может быть полезен в первую очередь как средство реконструкции моделей нейронов по экспериментальным данным внеклеточных или внутриклеточных записей от отдельных нейронов в мозге или в культуре. Полезность метода состоит в том, что рассматриваемая система демонстрирует большое богатство режимов, а предложенный подход ещё дополнительно снижает требования к точности модели за счёт табличной реконструкции нелинейной функции, то есть позволяет подогнать модель под более обширный класс экспериментальных сигналов. К тому же, метод работает по одной скалярной реализации переменной, соответствующей в модели трансмембранному потенциалу, то есть единственной реально доступной в эксперименте величине.
Библиографический список
References
Gouesbet G., Letellier C. Global vector-field reconstruction by using a multivariate polynomial L2 approximation on nets. Phys. Rev. E., 1994, vol. 49, pp. 4955-4972.
Smirnov D.A., Sysoev I.V., Seleznev E.P., Bezruchko B.P. Reconstructing nonautonomous system models with discrete spectrum of external action. Tech. Phys. Lett., 2003, vol. 29, no. 10, pp. 824-827.
Sysoev I.V., Prokhorov M.D., Ponomarenko V.I., Bezruchko B.P. Reconstruction of ensembles of coupled time-delay system from time series. Physical Review E, 2014, vol. 89, 062911. Belyustina, L.N., Shalfeev, V.D. Theory of a nonlinear system of automatic phase-frequency control. Radiophys Quantum Electron, 1968, vol. 11, pp. 213-220.
Mishchenko M.A. Neuron-like model based on a phase-locked loop. Journal of University of Nizhny Novgorod, 2011, vol. 5, no. 3, pp. 279-282 (in Russian).
Mishchenko M.A., Shalfeev V.D., Matrosov V.V. Neuron-like dynamics in phase-locked loop. Izvestiya VUZ. Applied Nonlinear Dynamics, 2012, vol. 20, no. 4, pp. 122-130. Matrosov V.V., Mishchenko M.A., Shalfeev V.D. Neuron-like dynamics of a phase-locked loop. The European Physical Journal. Special Topics, 2013, vol. 222, no. 10, pp. 2399-2405. Sysoev I.V., Sysoeva M.V., Ponomarenko V.I., Prokhorov M.D. Neuron-like dynamics in a phase-locked loop system with delayed feedback. Tech. Phys. Lett., 2020. In press. Kozlov A.K., Shalfeev V.D. Управление хаотическими колебаниями в генераторе с запаздывающей петлей фазовой автоподстройки. Izvestiya VUZ. Applied Nonlinear Dynamics, 1994, vol. 2, no. 2, pp. 36-48 (in Russian).
Bakunov G.M., Matrosov V.V., Shalfeev V.D. On quasi-synchronous regimes in a phase lock loop with the secondorder filter and approximate inclusion of the delay. Izvestiya VUZ. Applied Nonlinear Dynamics, 2011, vol. 19, no. 3, pp. 171-179 (in Russian).
Kashchenko S.A., Mayorov V.V. Modeli Volnovoy Pamyati. M.: Librocom, 2013, 288 p. (in Russian).
Glyzin S.D., Kolesov A.Y., Rozov N.K. Self-excited relaxation oscillations in networks of impulse neurons. Russian Mathematical Surveys, 2015, vol. 70, no. 3, pp. 383-452. Glyzin S.D., Kolesov A.Y., Rozov N.K. On a method for mathematical modeling of chemical synapses. Differential Equations, 2013, vol. 49, no. 10, pp. 1193-1210.
Glyzin S.D., Kolesov A.Y., Rozov N.K. Modeling the bursting effect in neuron systems. Mathematical Notes, 2013, vol. 93, no. 5, pp. 682-699.
Hodgkin A., Huxley A. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol., 1952, vol. 117, pp. 500—544. Mishchenko M.A., Bolshakov D.I., Matrosov V.V. Instrumental implementation of a neuronlike generator with spiking and bursting dynamics based on a phase-locked loop. Tech. Phys. Lett., 2017, vol. 43, pp. 596-599.
Sysoev I.V., Prokhorov M.D., Ponomarenko V.I., Bezruchko B.P. Determination of parameters of elements and coupling architecture in ensembles of coupled time-delay systems from their time series. Tech. Phys., 2014, vol. 59, no. 10, pp. 1434-1444.
Kapranov M.V. Capture band for phase-locked loop. Journal Radioengineering, 1956, vol. 11, no. 12, p. 37.
Belyustina L.N. The study of nonlinear phase-locked loop system. Radiophysics and Quantum Electronics, 1959, vol. 2, I. 2, pp. 277 (in Russian).
Sysoev I.V. Reconstruction of ensembles of generalized Van der Pol oscillators from vector time series. Physica D, 2018, vol. 384-385, pp. 1-11.
Сысоева Марина Вячеславовна - родилась в Саратове (1987). Окончила факультет нано- и биомедицинских технологий Саратовского государственного университета (2011). Защитила кандидатскую диссертацию по специальностям «Биофизика» и «Радиофизика» (2015). С 2015 года работает на кафедре «Радиоэлектроника и телекоммуникации» Саратовского государственного технического университета в должности доцента. Научные интересы -анализ временных рядов, нейронаука, математическое моделирование. Опубликовала свыше 20 научных статей по указанным направлениям.
Россия, 410019 Саратов, ул. Зелёная, 38
Саратовский филиал Института радиотехники и электроники
имени В.А. Котельникова РАН
Россия, 410054 Саратов, ул. Политехническая, 77
Саратовский государственный технический университет
имени Гагарина Ю.А.
E-mail: bobrichek@mail.ru
Сысоев Илья Вячеславович - родился в Саратове (1983). Окончил факультет нелинейных процессов СГУ (2004), защитил диссертацию на соискание учёной степени кандидата физико-математических наук (2007) и доктора физико-математических наук (2019). Профессор кафедры динамического моделирования и биомедицинской инженерии, ответственный секретарь редакционной коллегии журнала «Известия вузов. ПНД». Научные интересы - исследование сигналов биологической природы методами нелинейной динамики, исследование эффективности и модернизация подходов к анализу сигналов. Автор более 50 публикаций, индексируемых в международных базах.
Россия, 410019 Саратов, ул. Зелёная, 38
Саратовский филиал Института радиотехники и электроники имени В.А. Котельникова РАН Россия, 410012 Саратов, ул. Астраханская, 83 Саратовский национально исследовательский государственный университет имени Н.Г. Чернышевского Россия, 603950 Нижний Новгород, пр. Гагарина, 23 Национальный исследовательский Нижегородский государственный университет имени Н.И. Лобачевского E-mail: ivssci@gmail.com
Пономаренко Владимир Иванович - родился в Саратове (1960). Окончил Саратовский государственный университет (1982). Защитил диссертацию на соискание ученой степени кандидата физико-математических наук (1992) и доктора физико-математических наук (2008). Ведущий научный сотрудник Саратовского филиала Института радиотехники и электроники им. В.А. Котельникова РАН. Область научных интересов: нелинейная динамика, системы с запаздыванием, синхронизация, моделирование биологических систем. Автор более 200 научных публикаций.
Россия, 410019 Саратов, ул. Зелёная, 38
Саратовский филиал Института радиотехники и электроники имени В.А. Котельникова РАН Россия, 410012 Саратов, ул. Астраханская, 83 Саратовский национально исследовательский государственный университет имени Н.Г. Чернышевского E-mail: ponomarenkovi@gmail.com
Прохоров Михаил Дмитриевич - родился в Саратове (1968). Окончил Саратовский государственный университет (1992). Защитил диссертацию на соискание ученой степени кандидата физико-математических наук (1997) и доктора физико-математических наук (2008). Заведующий лабораторией Саратовского филиала Института радиотехники и электроники им. В.А. Котельникова РАН. Область научных интересов: нелинейная динамика и ее приложения, математическое моделирование, анализ временных рядов. Имеет более 200 научных публикаций.
Россия, 410019 Саратов, ул. Зелёная, 38
Саратовский филиал Института радиотехники и электроники имени В.А. Котельникова РАН E-mail: mdprokhorov@yandex.ru