Как мы создали наш сервис моделирования инфекционных заболеваний для Германии

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

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

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

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

Итак, без дальнейших церемоний, давайте начнем!

(Необязательно) Установка библиотеки CovPred

Если вы хотите поиграть сами, вы можете использовать библиотеку C++ из нашего репозитория GitHub:

GitHub — RefreshingAutoMate/EpidemicsRunnerAPI
API EpidemicsRunner — это библиотека C++, предоставляющая различные алгоритмы, полезные для моделирования инфекционных заболеваний…github.com

Библиотека EpidemicsRunnerAPI является преемницей библиотеки Epidemics, которая была плагином для программного обеспечения для моделирования. Однако Epidemics предлагает автономный графический интерфейс, который вы можете предпочесть. Он также был больше протестирован публикой, о чем можно судить по большему количеству звезд GitHub. Вы можете найти эту версию здесь:

GitHub — devanshr/Epidemics: подключаемый модуль Epidemics для программного обеспечения UG4
Epidemics — это подключаемый модуль UG4, предоставляющий различные алгоритмы для моделирования эпидемий. Этот документ должен дать широкий…github.com

Недавно мы представили привязки к Python 3. К моменту публикации этой статьи эти привязки могут стать общедоступными. Использование привязок Python — самый простой способ. Вы можете установить привязки Python через pip с помощью

$pip install covpred

В этой статье используются привязки Python.

Модель SEIRD

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

Тем не менее, эти симуляции по-прежнему полезны, поскольку они могут ответить на вопросы «что могло бы произойти, если бы наш местный регион удовлетворял аксиомам модели». Как мы выяснили, это часто так.

Модель (S)восприимчивого (E)выявленного (I)зараженного (R)обнаруженного (D)облегченного представляет собой разновидность почтенных компарментарных моделей, которые можно найти во многих учебниках по математике. Поскольку многие документы и статьи уже касаются теоретических описаний используемых алгоритмов, эта статья в основном сосредоточена на том, как можно использовать эти модели.

Этот конкретный вариант модели SEIRD был создан, в частности, Деваншем Растоги, одним из основателей нашего API EpidemicsRunner. В CovPred мы определяем классы как таковые:

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

Только люди из класса Exposed могут распространять болезнь в этой формуле. Предполагается, что люди из класса Infected остаются в полной изоляции.

Компартментальные модели этой формы не лишены критики, поскольку они предполагают, например, идеальное смешение всех классов S-E-I-R-D. На практике мы обнаружили, что их можно использовать, особенно в ограниченные сроки.

Как видно, эту модель можно настроить с помощью параметров. На скриншоте они называются alpha,kappa,pp,qq и theta. Мы находим следующее:

Альфа контролирует скорость перехода людей из класса «Восприимчивый» в класс «Выставленный».
Каппа контролирует процент людей, переходящих из класса «Подверженные» в класс «Зараженные и выздоровевшие».
qq контролирует скорость перехода людей из класса Exposed в класс Infected and Recovered.
pp контролирует скорость перехода людей из класса «Зараженные» в класс «Выздоровевшие и умершие».
тета контролирует процент людей, которые переходят из класса Зараженные в класс Выздоровевшие и Умершие

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

В CovPred мы определяем проблемы с начальными значениями для модели SEIRD ODE. Для этого теста в момент времени t=0 мы инициализируем нашу модель с помощью

S(0)=753056
E(0)=2714
I(0)=0
R(0)=0
D(0)=72

и симулировать на 44 дня.

В питоне это выглядит так:

import covpred
import matplotlib.pyplot as plt
import numpy as np
alpha=4.95954479066937e-07
kappa=0.356035567977659
theta=4.14932000304998e-07
qq=8
pp=5
u0=[753056,2714,0,0,72] #initial conditions
t_start=0 #simulation start time
t_end=43 #simulation end time
model=covpred.SEIRD(alpha,kappa,theta,qq,pp) #defining a SEIRD model
model.change_minimum_stepsize(0.01) #changing solver minimum stepsize
model.change_linear_implicit_maximum_stepsize(0.1)
timepoints, sim_data=model.run_linear_implicit(t_start,u0,t_end) #running the SEIRD model
sim_data=np.array(sim_data).reshape(len(timepoints),int(len(sim_data)/len(timepoints))) #reshape data so that result is not one long contiguous array anymore
plt.plot(timepoints,sim_data[:,0],label="Susceptibles")
plt.plot(timepoints,sim_data[:,1],label="Exposed")
plt.plot(timepoints,sim_data[:,2],label="Infected")
plt.plot(timepoints,sim_data[:,3],label="Recovered")
plt.plot(timepoints,sim_data[:,4],label="Deceased")
plt.legend()
plt.xlabel('time')
plt.ylabel('People')
plt.show()

Это дает:

Если теперь изменить альфу, то получатся следующие графики:

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

Как вы могли заметить в фрагменте кода, мы запустили модель с

timepoints, sim_data=model.run_linear_implicit(t_start,u0,t_end)

Вышеупомянутое означает, что система SEIRD ODE решается с помощью линейного неявного решателя. Этот решатель обладает хорошими свойствами стабильности, работает быстрее, чем полностью неявный решатель, и имеет автоматические свойства размера шага. Определяется здесь. Мы полагаем, что в MATLAB его эквивалент будет называться ODE23s.

Однако вы также можете запустить модель с помощью

timepoints, sim_data=model.run(t_start,u0,t_end)

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

В таких случаях, как на изображении выше, для обоих решателей требуется уменьшать размер шага до тех пор, пока не исчезнут дикие колебания. Мы обнаружили, что явный решатель быстрее приводит к осцилляциям и разрушению за конечное время. Тем не менее, оба решателя имеют свои места. Особенно, если требуется скорость (что мы позже увидим в формулировках PDE модели SEIRD), выбор между неявным и явным решателем может иметь большое значение.

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

Модель SEIRD с переменной альфой

На первый взгляд, дополнение модели SEIRD с помощью переменной альфы кажется очень простым:

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

Полную реализацию этой функции можно найти здесь. В EpidemicsRunner модели SEIRD с переменной альфа имеют суффикс _VARA. Давайте проверим это на Python!

import covpred
import matplotlib.pyplot as plt
import numpy as np
alphas=[1.95954479066937e-07,4.95954479066937e-7]
alpha_limits=[30] #change alpha values at this timepoint
kappa=0.356035567977659
theta=4.14932000304998e-07
qq=8
pp=5
u0=[753056,2714,0,0,72] #initial conditions
t_start=0 #simulation start time
t_end=43 #simulation end time
model=covpred.SEIRD_VARA(alphas, alpha_limits,kappa,theta,qq,pp)
model.change_minimum_stepsize(0.001)
model.change_linear_implicit_maximum_stepsize(0.01)
timepoints, sim_data=model.run_linear_implicit(t_start,u0,t_end) #running the SEIRD model
sim_data=np.array(sim_data).reshape(len(timepoints),int(len(sim_data)/len(timepoints))) #reshape data so that result is not one long contiguous array anymore
plt.plot(timepoints,sim_data[:,0],label="Susceptibles")
plt.plot(timepoints,sim_data[:,1],label="Exposed")
plt.plot(timepoints,sim_data[:,2],label="Infected")
plt.plot(timepoints,sim_data[:,3],label="Recovered")
plt.plot(timepoints,sim_data[:,4],label="Deceased")
plt.legend()
plt.xlabel('time')
plt.ylabel('People')
plt.show()

При t=30 значения альфа переключаются.

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

Модель SEIRD в ее текущей формулировке не выводит сводные данные для классов Susceptible, Exposed и Infected. Другими словами, графики Infected показывают, сколько людей одновременно инфицированы в момент времени t. В CovPred мы стараемся учитывать накопленные данные, а не параллельные заражения. Причина в том, что первое кажется проще. В свете наших уравнений это означает следующее:

Мы можем легко отслеживать наши кумулятивные незащищенные и зараженные, просто возвращая эти данные после запуска:

cumulated_exposed=model.get_cumulated_exposed()
cumulated_infected=model.get_cumulated_infected()

Для нашего предыдущего примера SEIRD_VARA это приводит к следующему графику:

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

Переменные альфа-значения также позволяют нам воспроизводить различные сценарии. Что может произойти, если через 10 дней будет обнаружен новый вариант, то есть альфа возрастет? Что может произойти, если будет введено строгое дистанцирование от контактов? Список возможных сценариев продолжается. В CovPred мы предоставляем пользователю различные такие сценарии.

В первой итерации службы CovPred мы использовали именно эту модель SEIRD, чтобы предоставить 429 различных оценок (1 страна, 16 Bundesländer и 412 Stadt-& Landkreise). Хотя результаты были обнадеживающими, это привело к одному большому основному предположению: а именно, что соседние регионы не взаимодействуют друг с другом.
Можно предположить, что это примерно верно для многих регионов. Тем не менее, у нас есть основания полагать, что такие регионы, как SK_Frankfurt_am_Main и SK_Offenbach, действительно влияют друг на друга.
Поэтому нам необходимо ввести некоторую форму связи между районами, чтобы отразить то обстоятельство, что некоторые регионы влияют друг на друга.

Связанная модель SEIRD

На основе произведения Деванша Растоги мы снова сели. В итоге мы остановились на следующей формулировке модели:

Это вводит две новые переменные, C, K и многие другие дополнительные параметры. Объяснения для классов прямолинейны.

C_i — это люди из класса «Восприимчивые» региона i, которые в настоящее время находятся в другом регионе.
K_i — это люди из незащищенного класса региона i, которые в настоящее время находятся в другом регионе.

Описания параметров также немного сложнее.

sigma_{i,j} влияет на количество восприимчивых и подвергшихся воздействию, которые перемещаются из области i в область j
omega_i влияет на количество восприимчивых лиц, покидающих область i
gamma_i влияет на скорость возвращения восприимчивых к региону i
phi_i влияет на скорость выставления покидающей области i
epsilon_i влияет на скорость возврата в регион i

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

Однако определение модели в Python намного проще. В следующем примере показано, как реализовать модель для двух связанных городов.

import covpred
import matplotlib.pyplot as plt
import numpy as np
class InitialConditions:
def __init__(self, u0, alpha,name):
  self.u0=u0
  self.alpha=alpha
  self.name=name
  self.kappa=0.5
  self.theta=0.1
  self.pp=1
  self.qq=1
  self.omega=2e-5
  self.gamma=0.1
  self.epsilon=1e-1
  self.phi=8e-3
alpha=[]
alpha_limits=[20]
kappa=[]
theta=[]
qq=[]
pp=[]
sigma=[]
omega=[]
gamma=[]
epsilon=[]
phi=[]
u0=[] #initial conditions of all regions
t_start=0
t_end=17
region_initials=[]
region_initials.append(InitialConditions([760000,0,0,0,0,0,0],[2.08e-6,0],"SK_Frankfurt_am_Main"))
region_initials.append(InitialConditions([160000,1500,0,0,0,0,0],[0,1e-6],"SK_Offenbach"))
for i in range(0,len(region_initials)):
  alpha+=region_initials[i].alpha
  kappa.append(region_initials[i].kappa)
  theta.append(region_initials[i].theta)        
  qq.append(region_initials[i].qq)
  pp.append(region_initials[i].pp)
  phi.append(region_initials[i].phi)
  epsilon.append(region_initials[i].epsilon)
  gamma.append(region_initials[i].gamma)
  u0+=region_initials[i].u0
  #Nobody from SK_Frankfurt_am_Main travels to SK_Offenbach
  if i==0:
    omega.append(0)
    for j in range(len(region_initials)):
      sigma.append(0)
  #But people from SK_Offenbach are traveling to SK_Frankfurt
  else:
    omega.append(1e-4)
    for j in range(len(region_initials)):
      sigma.append(1/(len(region_initials)-1))
model=covpred.SEIRD_CK(len(region_initials),alpha, alpha_limits,kappa,theta,qq,pp, sigma, omega, gamma,epsilon, phi)
model.change_minimum_stepsize(0.00001)
model.change_linear_implicit_maximum_stepsize(0.0001)
timepoints, sim_data=model.run_linear_implicit(t_start,u0,t_end)
sim_data=np.array(sim_data).reshape(len(timepoints),int(len(sim_data)/len(timepoints))) #reshape data so that result is not one long contiguous array anymore
fig, (ax1,ax2)=plt.subplots(2)
ax1.plot(timepoints,sim_data[:,0],label="Susceptibles")
ax1.plot(timepoints,sim_data[:,1],label="Exposed")
ax1.plot(timepoints,sim_data[:,2],label="Infected")
ax1.plot(timepoints,sim_data[:,3],label="Recovered")
ax1.plot(timepoints,sim_data[:,4],label="Deceased")
ax1.plot(timepoints,sim_data[:,6],label="Traveling Susceptibles from " +region_initials[0].name )
ax1.plot(timepoints,sim_data[:,7],label="Traveling Exposed from "+region_initials[0].name )
ax2.plot(timepoints,sim_data[:,5],label="Susceptibles")
ax2.plot(timepoints,sim_data[:,6],label="Exposed")
ax2.plot(timepoints,sim_data[:,7],label="Infected")
ax2.plot(timepoints,sim_data[:,8],label="Recovered")
ax2.plot(timepoints,sim_data[:,9],label="Deceased")
ax2.plot(timepoints,sim_data[:,10],label="Traveling Susceptibles from " +region_initials[1].name)
ax2.plot(timepoints,sim_data[:,11],label="Traveling Exposed from "+region_initials[1].name )
ax1.set_title(region_initials[0].name)
ax2.set_title(region_initials[1].name)
ax1.set_xlabel("time")
ax2.set_xlabel("time")
ax1.set_ylabel("People")
ax2.set_ylabel("People")
ax1.legend()
ax2.legend()
plt.show()

Как видно из начальных условий, в SK_Frankfurt_am_Main вначале не было зараженных или зараженных людей. Тем не менее, во Франкфурте происходят новые случаи из-за инфицированных людей, приезжающих из Оффенбаха.

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

Дальнейшие адаптации

Библиотека CovPred также включает модели SEIRD в формулировке дифференциальных уравнений в частных производных. Как показали Деванш Растоги и Марвин, эти модели в настоящее время не дают никаких преимуществ или даже хуже, чем наши связанные формулировки ОДУ, хотя для их решения требуется гораздо больше времени. Исследования этих моделей продолжаются с нашей стороны.

Будущее

Если вы заинтересованы в продолжении этой статьи, напишите нам письмо по адресу [email protected] и сообщите нам, что вы хотели бы услышать больше.

При наличии достаточного интереса мы покажем вам, как формулировать задачи оптимизации, решать их и анализировать результаты.

А пока берегите себя!