Прогнозирование уровня глюкозы в крови больных инсулинозависимым диабетом нейронными сетями и методом экстраполяции по выборке максимального подобия Текст научной статьи по специальности «Математика»

электронное научно-техническое издание

НАУКА и ОБРАЗОВАНИЕ

_Эл №ФС 77 - 30569. Государственная регистрация №0421000025. ISSH 1994-0408_

Прогнозирование уровня глюкозы в крови больных инсулинозависимым диабетом нейронными сетями и методом экстраполяции по выборке максимального подобия # 11, ноябрь 2010

авторы: Чернецов С. А., Чучуева И. А.

УДК.519.6

МГТУ им. Н.Э.Баумана serge.a. ch@gmail.com, chuchueva@yandex.ru

Введение

Сахарный диабет 1 типа - это метаболическое заболевание, вызванное абсолютным дефицитом секреции инсулина и характеризующееся неспособностью организма поддерживать уровень глюкозы в крови (BG - Blood Glucose) в целевом интервале 4-6 ммоль/л - в обычном состоянии и до 9 ммоль/л - после еды. Диабет вызывает множество опасных осложнений, избежать которые можно только путем контроля уровня BG и его удержания в физиологичном интервале. Основным путем решения этой задачи в настоящее время является введение в кровь пациента искусственных препаратов (генноинженерных человеческих инсулинов), которые могут симулировать действие эндогенного инсулина, вырабатываемого ß-клетками здоровой поджелудочной железы.

Оптимальные типы и дозы искусственного инсулина (далее - инсулина) зависят от многих факторов. Подбор этих типов и доз является сложной задачей, с которой могут справиться далеко не все пациенты. Для решения указанной задачи созданы системы непрерывного измерения уровня BG - Continuous Glucose Monitoring Systems (CGM-системы), а также системы непрерывного подкожного введения инсулина (инсулиновые помпы - insulin pumps).

На основе CGM-систем и инсулиновых помп разработаны и интенсивно разрабатываются системы автоматического управления уровнем глюкозы в крови пациента. С алгоритмической точки зрения эти системы включают в себя две следующие основные подсистемы: подсистема прогнозирования уровня BG; подсистема определения оптимального времени и требуемой дозы инсулина.

Известно значительное число публикаций, посвященных задаче прогнозирования уровня BG, среди которых можно выделить следующие три направления [1 - 3].

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

2) Использование нейросетевых алгоритмов прогнозирования значений уровня BG.

3) Использование комбинированных решений.

В рамках настоящей работы проводилось сравнение эффективности прогнозирования двух моделей - модель на основе нейронных сетей и модель экстраполяции временных рядов по выборке максимального подобия.

1. Постановка задачи

Пусть

Д = (Д^ Д 2 ) = (..., ^гр t_(n_2),..., ^ ^ tl, t2,..., ) -

дискретная временная сетка с постоянным шагом 5t = 20 мин. Здесь t0 - текущий момент времени; Д1 = (...^_(п_Х),t_(n_2),...,t_1,t0) - сетка «предыстории»; Д2 = t2,...,tm) - сетка, на которой выполняется прогноз.

Известны значения уровней BG и введенных пациенту доз инсулина в узлах сетки Д1 -векторы

Ь(Д1) = b_(n_2),..., b_1, Ь0)Т , /(Д1) = С., i_(n_1), i_(n_2),..., i_1, *0)Г .

Кроме того, известны значения компонентов вектора

с(Д1) = С - ^(п^ c_(n_2),..., c_1, С0)Т , где с_ 1 - количество углеводов принятых с пищей в момент времени t_j (оценивается и

задается пациентом).

Совокупность векторов Ъ(ДД/(ДД с(Дг) обозначим и(ДД

и = и(Д1) = (Ъ(ДД /^Д с(Д1)) = (...,ь_( П_1) ,..., Ъ0 ,... ,^'_(П_1) ,... ,..., c_(n_1),..., С0)Т .

Ставится задача прогноза уровня BG на период прогноза Т = [Хт -10], т.е. задача поиска

вектора у = Ъ(Д2) = (ЪьЪ2,...Ът)Т .

Критерии точности прогноза строятся на основе ошибок прогноза

=

] = 1,2,...,

где bj - фактическое значение уровня BG в момент времени tj. Точнее говоря, в качестве указанных критериев используются математическое ожидание mj ошибки прогноза Sj и ее среднеквадратичное отклонение оj.

Вектор фактических значений уровня BG в узлах сетки А 2 обозначим y = b(A2) = bb2,...)r .

2. Прогнозирование с помощью нейронных сетей

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

Было проведено исследование влияния на точность прогноза НС ряда параметров, которые можно разделить на две принцыпиальные группы - параметры сети и параметры исходного временного ряда:

• параметры сети (архитектура, алгоритм обучения, число нейронов скрытого слоя);

• параметры временного ряда (состав входных параметров, длина и характер линии задержки, длина - доля обучающей выборки).

Рассмотрены три типа нейронных сетей:

- Layer-Recurrent Network (далее LRN);

- сеть Элмана;

- Nonlinear Autoregressive Network (далее - NARX-сеть).

LRN имеют многослойную архитектуру, выход каждого из слоев (за исключением выходного) подается на его вход с временной задержкой на 1 шаг.

Сеть Элмана представляет собой частный случай LRN сети и имеет один внутренний слой, сигмоидальную функцию активации нейронов первого слоя, линейную функцию активации нейронов выходного слоя. Заметим, что MatLab R2009 поддерживает обучение сети Элмана только с помощью алгоритма градиентного спуска.

NARX-сеть является сетью прямого распространения, на вход сети подаются векторы

u,y, либо векторы u,y. Для формирования указанных векторов используются линии задержки TDL (Tapped Delay Line).

Для ускорения подготовки и проведения экспериментов разработан программный комплекс (рис. 1). Основу программы составляют классы dataHandler, Details, MatrixHost и NNmgr.

Класс dataHandler реализует однократное чтение данных из файла и создание матрицы данных - набора данных для различных интервала прогнозирования и количества учитываемых параметров.

Класс Details обеспечивает хранение информации об исследуемой сети или некотором наборе (матрице) таких сетей.

Рисунок 1. Структура MatLab-программного комплекса

Класс MatrixHost реализует обучение и тестирование указанного набора сетей. Объекты этого класса содержат матрицу объектов NNmgr, каждый из которых включает в себя подлежащую исследованию нейронную сеть, результат прогноза, а также ссылку на объект Details, который содержит все данные о матрице этих объектов. При создании указанной матрицы объектов задается список варьируемых параметров. Например, можно создать матрицу объектов, которая задает 4 различных число нейронов в скрытом слое, 5 алгоритмов обучения и 6 различных длительностей прогноза. Класс MatrixHost реализует также построение графиков результатов экспериментов.

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

3. Модель экстраполяции по выборке максимального подобия

Временной ряд Z=(z(t1), z(t2), z(t3),..., z(tN)) обозначим Zf =(z1, z2, z3,..., zN). Набор последовательных значений Zf =(zt, zt+i, zt+2,..., zt+M-i), лежащих внутри временного ряда

Z1W, назовем выборкой из этого ряда длины М с моментом начала отсчета М е [1: (Ы -1)]; ^ е[1:(N - М)]. Разность начал отсчетов выборок 2М , 21М_к будем называть задержкой k, k е [1: (7-1)].

Определим коэффициент корреляции выборок 2'М , как

Pk, (1)

уЩЁм]

где соу2м, ) - ковариация исходных выборок, а Ц[2М ] и ] - их дисперсии [4].

На основе коэффициента корреляции (1) введем меру подобия выборок 2'М , 2М

|соу2м , 2М_к )|

pM = £[0,1].

p[ZM ] ^D[Zf_k ]

Мера подобия рк зависит от длины M выборок Z t , Ztk, а также от задержки к, и характеризует степень линейной зависимости указанных выборок: чем ближе к единице значение ppf , тем выше их линейная зависимость.

Определим для выборки Zf величины р1, р2,..., p1f_1 и найдем их максимальное значение

M /MM M 1

Pkrnax=maX(Pi , Pi Pt-l). Выборку, соответствующую задержке kmax обозначим Z^max и назовем выборкой максимального подобия для выборки ZM .

Гипотеза подобия. Если выборки ZM , ZM,ii,r имеют значение величины pkiax близкое

a t J t—k max • k max

" r> Г rjM+P rjM+P M+P

к единице, то для некоторых значений P и выборок Zt , Zt—k max значение величины ркmax также близко к единице.

Перейдем к векторному обозначению выборки Zf =(zt, zt+1,..., zt+M-1)T и временного ряда

Z^ =(z1, z2,..., zN)T. Аппроксимируем выборку Zс помощью выборки ZM—kmax :

Z M = ai'Z Mk max +ao+E ;

Z M=arZ M—k max +ao. (2)

t—k г i M

Здесь а1, а0 - вещественные константы, Е - М-мерный вектор значений ошибок аппроксимации, а Ъ М аппроксимированные значения Z М .

Во введенных обозначениях выборку Z N+1 можно выразить через некоторую выборку Z р , лежащую внутри исходного временного ряда Z N , в виде

Z N+i = arZ P + ao. (3)

Алгоритм определения выборки Z р имеет следующий вид:

- для выборки Z M+1 находим выборку максимального подобия Z k^*, где

kmax *=N-M+ 1-kmax;

MM

- согласно формуле (2), аппроксимируем выборку Z N-M+1 с помощью выборки Z k max*:

Z M-M + 1 =ai'Z Mmax*+ao; (4)

P

- в соответствии с гипотезой подобия, в качестве выборки Z используем выборку

PP

ZT =Zkmax*+M, то есть выборку, расположенную на оси времени сразу за выборкой максимального подобия (рис. 2).

В формуле (4) вещественные константы a1, a0 определяются путем решения методом наименьших квадратов задачи

M

S = ^ (ZN-M+1-i - ZN-M+1-i) ^ min ,

i=1

которая детально рассмотрена в работе [5].

zM ZP 7M ZP

Zk max* Zk max*+M ZN-M+1 ZN+1

Рисунок 2. Положения выборок Zkfmax* , Zpmax*+M , ZN!-M+1 и ZL на оси времени

Таким образом, экстраполированные значения Z f+1 временного ряда Z f определяются формулой

Z f+1 =ai Z Pmax* +м +a0=EMMLS (M), (5)

которая представляет собой модель экстраполяции временных рядов по выборке максимального подобия (extrapolation model based on maximum likeness set, EMMLS). Разработанная модель имеет один параметр M. Методика идентификации модели детально рассмотрена в работах [6, 7].

4. Сравнение точности прогнозирования уровня BG с помощью нейронных сетей и

методом выборки максимального подобия

Рассматривался временной ряд уровней BG. Характеристики временного ряда приведены в таблице 1. Ставилась задача прогноза значений этого ряда на 20, 60 и 90 минут (на 4, 12, 18 значений временного ряда).

Для сравнения моделей был произведен прогноз контрольного периода длиной около 3800 значений на модели EMMLS (M) и модели искусственной нейронной сети (ИНС).

ИНС-модель представляла собой нейронную сеть прямого распространения, принимающая на вход пять параметров (ГК, инсулин, углеводы, гликемический индекс и физическую нагрузку) с использованием разреженно-суммирующей линии задержки длинной, раной три (т.е. общее число входов равно 15). Сеть обучалась алгоритмом Левенберга-Маркадта.

Таблица 1. Параметры используемого временного ряда

Шаг по времени Временной ряд Длина ряда Среднее значение Стандартное отклонение Минимум Максимум

300 с. BG, ммоль/л 29 640 2.25 14.24 2.20 21.90

Сравнение результатов прогнозирования моделей приведено в таблице 2.

Краткосрочный прогноз на моделях EMMLS (18), EMMLS (78) для времени упреждения двадцать минут и час, соответственно, в среднем имеет точность выше аналогичного прогноза с использованием ИНС. Как видно из таблицы 2, значение МАРЕ для модели EMMLS (18) при Р=4 составляет 5,07%, а при Р=12 - 11,33%. Аналогично, значение

МАРЕ для модели ИНС при Р=4 составляют 8,09%, а при Р=12 - 13,21%. Важно отметить, что результаты прогнозирования позволяют выявить сильные изменения уровня BG (резкие увеличения и снижения этого уровня).

Таблица 2. Точность прогноза моделей ИНС и EMMLS (М)

Модель Р MAE, ммоль/л MAPE, % Число точек точнее, % Период экстраполяции, час Период идентификации, час

EMMLS (18) 4 0.36 5.07 45 0.35 0.5

ИНС 0.30 4.12 55 80 80

EMMLS (78) 12 0.79 11.33 52 0.52 0.85

ИНС 0.91 12.21 48 80 80

EMMLS (180) 18 0.97 14.70 58 0.10 2.50

ИНС 1.12 17.02 42 80 80

При долгосрочном прогнозировании (на полтора часа) модель EMMLS (М) показала точность более, чем на 2% превышающую точность, достигнутую с помощью ИНС.

Была также исследована эффективность совместного использования ИНС и EMMLS(M) для долгосрочного прогнозирования. Рассмотрена линейная комбинация вида

V 1 V 1 V 1

Zagg N+1 =0.5416^ем1^8 N+1 +0.5808^АЖ N+1 -0.7326,

11

где Ъемм^ +1 - прогноз на модели EMMLS(M), +1 - прогноз с использованием ИНС,

1

а Zagg +1 - итоговое значение прогноза BG. Исследование показало, что такая комбинация обеспечивает значение MAPE, равное 12.13 %, а значение MAE, равное 0.77 ммоль/л.

Заключение

В статье выполнено сравнение эффективности прогнозирования уровня BG с использованием нейронных сетей и модели EMMLS. Показано, что при краткосрочном прогнозировании уровня BG более точный результат дает нейронная сеть, а при долгосрочном прогнозировании - модель EMMLS.

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

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

Работа выполнена при поддержке гранта РФФИ №10-08-00816-а.

Литература

1. Гоменюк С.М., Емельянов А.О., Карпенко А.П., Чернецов С.А. Методы прогнозирования оптимальных доз инсулина для больных сахарным диабетом I типа. Обзор // Наука и образование: электронное научно- техническое издание, www.technomag.edu.ru, апрель, 2009.

2. Гоменюк С.М., Емельянов А.О., Карпенко А.П., Чернецов С.А. Обзор методов и систем прогнозирования оптимальных доз инсулина для больных сахарным диабетом I типа // Информационные технологии. 2010 (в печати).

3. Chernetsov S.A., Emelyanov A.O., Karpenko A,P. Research of Neural Network-Based Blood Glucose Level Forecasting Systems for Insulin-Dependant Diabetes Patients // The 6th international workshop on Wearable Micro and Nanosystems for Personalised Health (in print).

4. Бокс Дж., Дженкинс Г.М. Анализ временных рядов, прогноз и управление. - М.: Мир, 1974. - С. 406.

5. Draper N. R., Smith H. Applied regression analysis. - New York: Wiley, In press, 1981. - С. 693.

6. Павлов Ю. Н., Чучуева И. А. Экстраполяция псевдослучайных процессов по максимуму подобия. [Электронный ресурс] // Электронное научно-техническое издание «Наука и образование». - 2009 г - № 7. (http://technomag.edu.ru/doc/129712.html)

7. Чучуева И.А. Модель экстраполяции по максимуму подобия (ЭМП) для временных рядов цен и объемов на рынке на сутки вперед ОРЭМ (Оптовом рынке электроэнергии и мощности). [Электронный ресурс] // Электронное научно-техническое издание «Наука и образование». - 2010 г - № 1. (http://technomag.edu.ru/doc/135870.html)