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

I. Получение минимальной модели

Фармакокинетика изучает всасывание, распределение и метаболизм лекарственных средств в организме. Минимальная модель — это фармакокинетическая модель, предложенная Bergman et al. в 1979 году и с тех пор был изменен и расширен.²

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

Бергман и его соавторы попробовали несколько систем для моделирования взаимодействия глюкозы и инсулина, прежде чем остановились на окончательной версии модели. Они попытались разработать «самую простую модель, основанную на известной физиологии, которая могла бы объяснить связь между инсулином и глюкозой, обнаруженную в данных» (Bergman et al, 1979), и опробовали семь вариантов модели.

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

  • G - концентрация глюкозы в крови как функция времени
  • Х - концентрация инсулина в тканевой жидкости в зависимости от времени.
  • I - концентрация инсулина в крови, измеренная при вводе
  • Gᵦ - равновесная (базальная) концентрация глюкозы в крови, измеренная до исследования
  • Iᵦ - равновесная (базальная) концентрация инсулина в крови, измеренная до исследования
  • pᵢ - параметры, контролирующие скорость создания/удаления глюкозы и инсулина

Гормональная система организма пытается поддерживать гомеостаз таким образом, чтобы концентрация глюкозы в крови поддерживалась на равновесном уровне. Минимальная модель отражает тот факт, что изменения уровня глюкозы dG/dt и уровня инсулина в тканевой жидкости dX/dt стремятся вернуть организм к гомеостазу.

Гипергликемия возникает, когда поджелудочная железа не вырабатывает достаточного количества инсулина (диабет I типа) или если рецепторы инсулина становятся нечувствительными (диабет II типа). Тяжелая гипергликемия является ключевым симптомом сахарного диабета, распространенного и серьезного заболевания, поражающего значительную часть населения США². Минимальная модель дает представление о функции поджелудочной железы человека через реакцию инсулина на инъекции глюкозы.

II. Физико-информированные нейронные сети

Нейронная сеть, основанная на физике (PINN) — это аппроксиматор функций, который объясняет физические законы Вселенной с помощью языка дифференциальных уравнений. Эти типы сетей процветают в средах с низкими ограничениями данных, в ситуациях, когда типичные подходы машинного обучения дают сбои. Используя дифференциальные уравнения для управления обучением нейронной сети, мы можем получить приемлемое решение для нашей системы, несмотря на отсутствие существенных данных.

Существует множество решателей дифференциальных уравнений, которые хорошо подходят для аппроксимации решений динамических систем. Они часто дают быстрые решения для многих систем дифференциальных уравнений. Так зачем использовать нейронную сеть, когда доступны обычные решатели дифференциальных уравнений? Есть две области, в которых PINN выгодны:

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

Как работает PINN?

Во-первых, давайте определим нейронную сеть как функцию u(x; θ), которая отображает вход x в выход y, где θ представляет веса и смещения сети. Цель обучения нейронной сети — найти оптимальный набор весов и смещений θ*, который минимизирует ошибку прогнозирования между выходными данными модели y и истинными выходными даннымиy*.

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

Пусть u₁(t; θ₁) представляет решение нейронной сети для изменения уровня глюкозы относительно времени (dX/dt) и u₂ (t; θ₂) представляют решение нейронной сети для изменения концентрации инсулина в тканевой жидкости относительно время (dX/dt). θ₁ и θ₂ представляют наборы весов и смещений для каждой соответствующей функции решения.

Цель состоит в том, чтобы найти θ₁* и θ₂* такие, что невязка между производной выхода функции решения, назовем ее uₙₙ и соответствующая вынуждающая функция (dG/dt, dX/dt) минимизируется. Идея, лежащая в основе этого, заключается в том, что, если производная выхода нейронной сети соответствует дифференциальному уравнению, решение нейронной сети будет подчиняться законам динамической системы.

Определим дифференциальный оператор L. Мы можем взять производные от выходной функции решения, u₁,ₙₙи u₂,ₙₙ, используя L, поскольку нейронная сеть состоит из дифференцируемых функций активации. Тогда остатки определяются следующим образом:

Затем это можно выразить как задачу оптимизации с использованием L² потери остатка:

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

(Примечание: переключение между числовыми нижними индексами 1, 2 и нижними индексами G, X было вызвано тем, что на medium.com нет возможности набирать нижние индексы без их копирования и вставки. Некоторые буквы по какой-то причине не могут быть скопированы.)

III. Решение минимальной модели

Давайте попробуем эту технику на практике. Мы будем использовать данные о концентрации глюкозы/инсулина, полученные от Pacini et al. (1998)⁴, который затем был размещен в сети благодаря Аллену Дауни²:

from os.path import basename, exists

def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve
        local, _ = urlretrieve(url, filename)
        print('Downloaded ' + local)

download('https://github.com/AllenDowney/ModSim/raw/main/data/' +
         'glucose_insulin.csv')

'https://github.com/AllenDowney/ModSim/raw/main/data/glucose_insulin.csv'

Данные, которые мы рассматриваем, показывают концентрацию глюкозы и инсулина в крови одного пациента с течением времени. Выполняемый тест называется внутривенным тестом на толерантность к глюкозе с частым отбором проб (FSIGTT), при котором субъекту натощак вводят глюкозу, а образцы крови собирают с интервалами 2–10 минут в течение 3 часов.

Нам нужно импортировать необходимые пакеты для начального анализа и считать данные. Библиотека для настройки PINN будет представлена ​​позже.

import numpy as np
import pandas as pd

from scipy.interpolate import interp1d
from scipy.interpolate import UnivariateSpline
from scipy.integrate import solve_ivp
from scipy.optimize import leastsq

import seaborn as sns
import matplotlib.pyplot as plt
df = pd.read_csv('glucose_insulin.csv', index_col='time')
df.head()

Для краткости некоторые детали будут опущены, на которые можно ссылаться в полном блокноте ниже. У нас есть дискретный набор точек данных, и нам нужно сделать данные гладкими и непрерывными, чтобы мы могли позже производить выборку из точек. Мы можем использовать функцию одномерного разделения scipy для выполнения этой задачи:

t_0 = df.index[0]
t_f = df.index[-1]
t = np.linspace(t_0, t_f, int(t_f-t_0))

G_func = UnivariateSpline(df.glucose.index[1:], df.glucose.values[1:], k=2, s=100)
G_series = pd.Series(G_func(t), index=t)

I = UnivariateSpline(df.insulin.index, df.insulin.values, k=1, s=200)
I_series = pd.Series(I(t), index=t)

Теперь мы можем построить данные, чтобы увидеть, как измерялись уровни пациентов:

На высоком уровне мы можем видеть, что после инъекции глюкозы концентрация инсулина в крови I(t)имеет слегка запаздывающий всплеск от глюкозы в ответ на инъекцию. По мере того, как инсулин выполняет свою работу, уровни глюкозы и инсулина снижаются. Инсулин будет диффундировать между кровью и тканевой жидкостью по мере необходимости, постепенно снижаясь. Также действует естественное снижение уровня глюкозы.

Нам нужно определить минимальную модель в коде. Прежде чем решать систему ОДУ с помощью PINN, давайте попробуем подход с численными методами. Мы будем использовать функцию задачи Scipy с начальным значением, которая использует метод аппроксимации под названием Рунге-Кутта-Дорманд-Принс (RKDP).

def minimal_model(t, states, p1, p2, p3):

    G = states[0]
    X = states[1]

    dGdt = -p1 * (G - Gb) - X*G
    dXdt = p3 * (I(t) - Ib) - p2 * X

    return [dGdt, dXdt]

Теперь, когда мы определили нашу модель, нам нужно найти параметры, которые соответствуют данным функции. Мы можем определить функцию ошибки и использовать метод наименьших квадратов, чтобы минимизировать потери между выходом RKDP и данными. Результат обеспечит наши оптимальные параметры.

def error(parameters, df):
    G0, X0 = 270, 0
    p1, p2, p3 = parameters
    sol = solve_ivp(minimal_model, t_span=(t_0,t_f), y0=[G0, X0], args=(p1,p2,p3), method='RK45', dense_output=True)
    errors = sol.sol(df.glucose.index)[0]- df.glucose.values
    # return errors starting at index 2 to allow insulin to catch up to glucose injection
    return errors[3:]

# start with initial value estimates
param_estimates = [0.01, 0.01, 1.0e-05]
# use scipy's leastsq function to get optimal parameter values (takes a while
# to run)
optimal_params = leastsq(error, x0=param_estimates, args=(df))[0]
print(optimal_params)

Найдя оптимальные параметры, мы можем использовать метод RKDP для решения системы и построить результаты:

mm_glucose, mm_insulin = solve_ivp(minimal_model, t_span=(t_0,t_f), y0=[G0, X0], args=optimal_params, method='RK45', dense_output=True)

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

Теперь давайте попробуем использовать PINN для моделирования той же системы. Мы используем пакет под названием neurodiffeq, который поможет с настройкой нейросетей и оптимизаторов. Это необходимые функции:

from neurodiffeq import diff
from neurodiffeq.solvers import Solver1D
from neurodiffeq.monitors import Monitor1D
from neurodiffeq.conditions import IVP
from neurodiffeq.networks import FCNN, SinActv
from neurodiffeq.ode import solve_system
from neurodiffeq.generators import Generator1D

import torch
import torch.nn as nn

Затем мы снова настраиваем минимальную модель, но на этот раз мы включаем дифференциальный оператор, который будет составлять нашу функцию потерь. В этом подходе нейронные сети имеют два скрытых слоя по 32 единицы с активационными функциями. Данные обучения и проверки будут единообразно отбираться из нашей функции интерполяции ранее, и мы прогоняем 20000 эпох. Это займет несколько минут.

monitor = Monitor1D(t_min=t_0, t_max=t_f, check_every=1000)
train_gen = Generator1D(size=64,  t_min=t_0, t_max=t_f, method='uniform')
valid_gen = Generator1D(size=64, t_min=t_0, t_max=t_f, method='uniform')

def minimal_model_nn(G, X, tspan):
    dGdt = -p1*(G - Gb) - X*G
    dXdt = (torch.tensor(p3*(I(tspan.detach().numpy()) - Ib)) - p2*X)*10**4

    return [diff(G,tspan) - dGdt, diff(X,tspan) - dXdt]

conditions = [IVP(t_0=t_0, u_0=G0), IVP(t_0=t_0, u_0=X0)]
nets = [FCNN(n_input_units=1, n_output_units=1, hidden_units=(32, 32), actv=nn.Tanh),\
        FCNN(n_input_units=1, n_output_units=1, hidden_units=(32, 32), actv=nn.Tanh)]

solver = Solver1D(minimal_model_nn
                  , conditions
                  , t_min=t_0
                  , t_max=t_f
                  , nets=nets
                  , train_generator=train_gen
                  , valid_generator=valid_gen)

solver.fit(max_epochs=20000, monitor=monitor)
solution = solver.get_solution()

# The solution variable stores the trained neural net functions
# To plot the functions, we can evaluate them over our given time interval
G_nn, X_nn = solution(t, to_numpy=True)

Теперь мы можем построить функции решений:

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

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

IV. Обсуждение

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

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

Гилпин и др. (2020) предлагают тандемный метод моделирования систем с большими данными. Подход к обучению без учителя используется для извлечения динамики из больших наборов данных. В свою очередь, извлеченная динамика используется для разработки механистических математических моделей. Модель машинного обучения можно анализировать с помощью инструментов из теории динамических систем, получая преимущества обоих подходов³. Подход PINN представляет собой своего рода гибрид моделирования сверху вниз и снизу вверх. Динамическая система получается путем анализа динамики отдельных лиц в системе (снизу вверх), а затем эта динамика соответствует эмпирической модели (сверху вниз).

В данном случае PINN оказался не таким успешным, как численный метод. Есть несколько причин, которые кажутся правдоподобными. Для начала мы работаем с зашумленными реальными данными, собранными вручную от человека. Несмотря на наши усилия по сглаживанию данных при интерполяции, нейронная сеть производит выборку из данных с большим выбросом, что может привести к тому, что нейронная сеть узнает массивный фантомный выброс в данных. Нейронная сеть, которую мы пробовали, представляла собой стандартную готовую архитектуру. Возможно, с более тонкой настройкой наши результаты выглядели бы более заманчиво. Более интересными приложениями PINN могут быть системы, в которых численные решатели дают сбои, возможно, в более сложных многомерных системах.

При этом я надеюсь, что читатель сможет (1) лучше понять PINN, (2) как реализовать его на практике и (3) как выглядит моделирование реальных несовершенных данных.

V. Ссылки

[1]. Бергман, Р. Н. и др. (1979). Количественная оценка чувствительности к инсулину. Американский журнал физиологии vol. 236,6 стр. 667–77. doi:10.1152/ajpendo.1979.236.6.E667.

[2]. Дауни, Аллен (2022). Фармакокинетика. Моделирование и имитация в Python, https://allendowney.github.io/ModSimPy/chap17.html.

[3]. Гилпин, Уильям и др. (2020). Динамика обучения на основе больших наборов биологических данных: машинное обучение и системная биология. Современное мнение в области системной биологии, Vol. 22, с. 1–7, ISSN 2452–3100, https://doi.org/10.1016/j.coisb.2020.07.009.

[4]. Пачини, Джованни и др. (1998). Чувствительность к инсулину и эффективность глюкозы: минимальный модельный анализ обычного и модифицированного инсулином FSIGT. Американский журнал физиологии. 274. Е592–9. doi: 10.1152/ajpendo.1998.274.4.E592.