Сетевые модели эпидемии

Введение

Самая первая публикация, посвященная математическому моделированию эпидемий, датируется 1766 годом. В этой основополагающей статье Essai d'une nouvelle, посвященной анализу смертности среди детей, Даниэль Бернулли разработал математическую модель для анализа смертности от оспы. В дальнейшем, Хамер и Росс, применив одними из первых математическое описание распространения заболеваний, смогли решить задачи по выяснению механизмов регулярного повторения эпидемии кори и установлению взаимосвязи между количеством комаров и возникновением малярии. Труды Хамера и Росса, наряду с более поздними работами Росса и Хадсона, Сопера, а также Кермака и Маккендрика (Kermack, McKendrick, 1927), послужили теоретической основой для дальнейших исследований в области математического моделирования эпидемий.
В упомянутых работах для задач эпидемиологии был впервые применён так называемый «закон действующих масс», согласно которому количество вновь инфицированных в популяции прямо пропорционально произведению текущих численностей восприимчивых и инфицированных индивидов. Модель Кермака и Маккендрика дала начало широкому применению детерминированных SIR–моделей («Susceptible — Infected — Recovered» : Восприимчивый – Зараженный – Выздоровевший), в которых с помощью систем дифференциальных уравнений (непрерывное время) или разностных уравнений (дискретное время) описывается динамика групп восприимчивых, инфицированных и выздоровевших индивидов.
В 1957 году Кендаллом была сформулирована одна из первых пространственных моделей распространения эпидемий на основе уравнений в частных производных (Kendall, 1957). В том же году Бартлетт описывает распространение эпидемии на узлах пространственной структуры размера 6×6, пользуясь методами имитационного моделирования — нового направления, основанного на компьютерном моделировании изучаемых систем с применением численных методов Монте–Карло. В 1971 году Фоксом и Элвбэк с соавторами была опубликована первая индивидуум–ориентированная модель распространения заболевания (Elveback, 1971).
В некоторых задачах математической эпидемиологии требуется учитывать географическую и социальную неоднородность распространения инфекции. В этом случае используются имитационные модели, в которые включаются пространственные структуры. Ими могут быть клеточные автоматы и сетевые модели. В этих случаях рассматривается не пространственная близость индивидов популяции, а их близость в социальном отношении (социальные связи).


Базовая модель динамики взаимодействия иммунной системы человека с вирусом
Биология процесса и постановка задачи

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

5391141943735Неинфицированные клетки (х)
0Неинфицированные клетки (х)
24631652034540Вирус
0Вирус
33966151939290Инфицированные клетки (у)
0Инфицированные клетки (у)


Рисунок 1. Базовая модель вирусной динамики. Неинфицированные клетки «реагируют» со свободным вирусом, вызывая заражение клеток со скоростью β. Зараженные клетки производят свободные вирионы со скоростью k. Неинфицированные клетки, свободный вирус и инфицированные клетки погибают со скоростью d, u и a соответственно. Неинфицированные клетки пополняются с постоянной скоростью λ.

Свободные вирусные частицы заражают неинфицированные клетки со скоростью, пропорциональной произведению их численности, βxv. Константа скорости, β, описывает эффективность этого процесса, включая скорость, с которой вирусные частицы обнаруживают незараженные клетки, скорость проникновения вируса, а также скорость и вероятность успешного заражения.
Зараженные клетки продуцируют свободный вирус со скоростью, пропорциональной их численности, ky. Зараженные клетки погибают со скоростью ay, а свободные вирусные частицы удаляются из системы со скоростью uv. Следовательно, средняя продолжительность жизни инфицированной клетки составляет 1 / а, тогда как среднее время жизни свободной вирусной частицы составляет 1 / u. Общее количество вирусных частиц, образующихся в одной зараженной клетке, «размер разрыва», составляет к / а.
Неинфицированные клетки продуцируются с постоянной скоростью λ и погибают со скоростью dx. Среднее время жизни неинфицированной клетки составляет 1 / d. При отсутствии инфекции динамика популяции клеток-хозяев определяется как ẋ (t) = λ - dx. Это простое линейное дифференциальное уравнение. Без вируса количество незараженных клеток сходится к равновесному значению λ / d.
Объединяя динамику вирусной инфекции и клетки-хозяина, получаем базовую модель вирусной динамики (1), которая описывает вирусную динамику ВИЧ-1, и определяется следующим набором дифференциальных уравнений:
ẋ (t) = λ - dx - βxv,
ẏ(t) = βxv - ay, (1)
v̇ (t) = ky - uv.
Здесь x представляет количество неинфицированных клеток;
y количество инфицированных клеток;
v количество вирусных частиц в момент времени t;
λ - общая скорость производства новых неинфицированных клеток в единицу времени;
d, a и u - показатели смертности неинфицированных клеток, инфицированных клеток и вирусных частиц соответственно.
Это система нелинейных дифференциальных уравнений.
Аналитическое решение временных переменных не представляется возможным, но возможно вывести различные приближения и тем самым получить более полное понимание системы.
До заражения имеем y=0, v=0, а неинфицированные клетки находятся в равновесии x= λ / d. Обозначим через t=0 время возникновения инфекции. Предположим, что заражение происходит с определенным количеством вирусных частиц, v0. Таким образом, начальные условия x= λ / d, y=0 и v0. Способность вируса расти и вызывать инфекцию зависит от состояния, очень похожего на распространение инфекционного заболевания в популяции множества людей.
Решающее значение имеет базовое репродуктивное отношение R0, которое определяется как число вновь инфицированных клеток, возникающих из любой инфицированной клетки, когда почти все клетки не заражены.
Каждая зараженная клетка производит свободные вирусные частицы со скоростью k в единицу времени. Зараженные клетки получают из неинфицированных клеток и свободных вирусных частиц с общей скоростью βxv в единицу времени. Базовое репродуктивное число
R0 = βλk adu
как ожидаемое число вторичных вирусных частиц, продуцируемых одной новой вирусной частицей, попадающей в популяцию неинфицированных клеток в равновесии, в которой нет ни инфицированных клеток, ни вируса.
Если R0 <1, то вирус не будет распространяться, поскольку каждая зараженная клетка в среднем будет производить меньше, чем одна другая зараженная клетка. Цепная реакция является докритической. В среднем ожидается 1 / (1- R0) раундов заражения, прежде чем вирусная популяция вымирает.
Если, с другой стороны, R0> 1, то каждая зараженная клетка в среднем будет производить более одной вновь зараженной клетки. Цепная реакция вызовет взрывное размножение вируса. Рост вируса не будет продолжаться бесконечно, поскольку количество незараженных клеток ограничено. Там будет пик вирусной нагрузки и впоследствии затухание колебаний до равновесия.
Для этой модели существуют два равновесия. Равновесие без болезней (DFE) равно
x = x̂ = λ /d
y = ŷ = 0,
v = v̂ = 0
и эндемическое равновесие
x = x* = au βk,
y = y* = (R0 - 1) duβk,
v = v* = (R0 - 1) dβ.
В состоянии равновесия любая зараженная клетка в среднем приводит к появлению одной вновь зараженной клетки

. Таким образом, доля свободных вирусных частиц, которым удается заразить новые клетки, дается обратной величиной, а/k. Вероятность того, что клетка (рожденная незараженной) останется незараженной в течение ее жизненного цикла, составляет 1 / R0. Соотношение равновесия неинфицированных клеток до и после заражения составляет x0 / x* =R0.
Если вирус имеет базовое репродуктивное отношение намного больше единицы, тогда x* будет значительно уменьшено по сравнению с x0, что означает, что во время заражения равновесное количество неинфицированных клеток намного меньше, чем до заражения. Другими словами, приведенная выше простая модель не может объяснить ситуацию, когда во время персистирующей вирусной инфекции почти все заражаемые клетки остаются незараженными (x *≈x0), за исключением случая, когда R0 лишь немного больше единицы (что априори маловероятно в большинстве случаев).
Кроме того, если R0≫ 1, то равновесное количество инфицированных клеток и свободного вируса приблизительно определяется как y* ≈ λ / а и v*≈ (λk / au). Интересно, что обе величины не зависят от параметра инфекции β. Причина в том, что высокоинфекционный вирус (большее β) быстро заразит неинфицированные клетки, но в равновесии в системе будет только несколько неинфицированных клеток. Менее инфекционный вирус (меньшее β) займет больше времени, чтобы заразить неинфицированные клетки, но равновесное количество незараженных клеток выше. Для обоих вирусов произведение βx будет одинаковым в равновесном состоянии, что приведет к постоянной скорости продукции новых инфицированных клеток и, следовательно, к сходному равновесному содержанию инфицированных клеток и свободного вируса.
Для высоко цитопатического вируса (a намного больше, чем d) равновесное количество инфицированных клеток будет небольшим по сравнению с количеством клеток до заражения. Фактически, чем больше a, тем меньше численность как инфицированных клеток, так и свободного вируса. Для нецитопатического вируса (a~d) равновесное количество инфицированных клеток будет примерно эквивалентно общему количеству восприимчивых клеток до заражения.

Моделирование вирусной динамики

Рассмотрим развитие эпидемии, сходной с H1N1. Население будет разделено на три отдельные группы восприимчивых, зараженных и выздоровевших. Хотелось бы точно определить, какие параметры будут контролировать или предотвращать эпидемию. Население будет моделироваться тремя дифференциальными уравнениями. Команду Matlab ode45 можно использовать для решения таких систем дифференциальных уравнений.
Сделаем следующие основные допущения.
Допустим, что население представляет собой непересекающееся множество
Восприимчивых (susceptible) = S (t),
Зараженных ( Infected) = I (t) ,
Выздоровевших ( Recovered) = R (t).
Тогда общая численность населения N = S (t) + I (t) + R (t).
Допустим, все зараженные в конце концов выздоровеют.
Допустим, все выздоровевшие более не подвержены инфекции.
Допустим, что увеличение числа инфицированных пропорционально изменению во времени числа зараженных и числа восприимчивых.
Изменение в зараженной популяции будет увеличиваться из восприимчивой группы и уменьшаться в восстановленной группе.
I (t + Δt) – I (t) = Δt α S (t) I (t) – Δt β I (t),
где α отражает степень заразности или инфекционности эпидемии и β отражает скорость выздоровления. Теперь разделим на Δt и устремив его к нулю, чтобы получить дифференциальное уравнение для I:
dIdt= α S I - β I, I0=i0.
Дифференциальные уравнения для S и R получаются аналогичным образом:
dSdt=- α S I, S0=s0 .
dRdt= β I, r0=0.
Тогда скрипт модели эпидемии SIR:
% ввод времени и начальных условий и получение выходного вектора
% dS/dt, dI/dt, dR/dt
%Function sir = sir(t, x)
% Переменные
% альфа = степень заразности
alpha = .9;
% beta = скорость выздоровления
% tpop = общее количество населения в субпопуляции i
tpop = 11;
% ds/dt = sir(1)
% x(1) = восприимчивое население
ds = - (alpha*x(1)*x(2)/tpop);
% di / dt = sir(2)
% x(2) = зараженные
di = (alpha*x(1)*x(2)/tpop) - (beta*x(2));
% dr/ dt = sir(3)
% x(3) = выздоровевшие
dr = (beta*x(2));
sir = [ds, di, dr];
%
function sir = sir(t,x)
%Variables
%alpha=the global transmission rate
alpha=.3;
%beta=contact rate per day by people in subpopulation i
beta=.9;
%gamma=Infectious mortality of subpopulation i, from infected class
%gamma=.01;
%tow=Rate of progression from exposed to infected (? ?1 is latent period )
tow=1;
%row=Recovery rate
row=.5;
%tpop=Total number of people in subpopulation i
tpop=10;
%ds/dt=sir(1)
%x(1)=suscepatiblepeople
ds=-(beta*alpha*(((x(2)+x(3))*x(1))/tpop));
%de/dt=sir(2)
%x(2)=exposed by H1N1
%di=(beta*alpha*(((x(2)+x(3))*x(1))/tpop))-(tow*x(2));
%di/dt=sir(3)
%x(3)=Infected by H1N1
di=(alpha*(((x(2)+x(3))*x(1))/tpop))-(tow*x(2));
%dr/dt=sir(4)
%x(4)=Recovered from H1N1
dr=(row*x(2));

sir=[ds;di;dr];


Команда Matlab ode45 используется для приближенного решения системы дифференциальных уравнений. Эта команда является надежной реализацией для систем дифференциальных уравнений, в которой используется метод переменного размера шага и метод Рунге-Кутты четвертого и пятого порядка.
Поскольку существует три неизвестных и три дифференциальных уравнения, и нужно использовать функцию ode45 Matlab, файл вызова должен содержать три правые части дифференциальных уравнений, где y (1) = S, y (2) = I и y (3) = R. Выбор a и b в sir.m будет определять, будет ли I (t) увеличиваться; если I '(0) = (a S (0) - b) I (0)> 0, тогда I (т) увеличится. Файл sir.m содержит вызов ode45, и графики трех групп населения генерируются с помощью командного графика Matlab. Начальные популяции S (0) = 99, I (0) = 1 и R (0) = 0.
%
[t, x] = ode45(@sir, [0, 5], [100, 1, 0]);
p1 = plot(t, x);
xlabel('время');
ylabel('население'), legend('восприимчивые', 'зараженные', 'выздоровевшие');
>> grid on
[t, x] = ode45(@sir2, [0, 10], [100, 1, 0]);
p1 = plot(t, x);
xlabel('время');
ylabel('население'), legend('восприимчивые', 'зараженные', 'выздоровевшие');
>> grid on
Результатом работы программы является график:

Рисунок 2. Динамика распространения инфекционного заболевания

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

Детерминированная модель распространения ВИЧ

Аналитические результаты предполагают существование двух равновесий для системы (1) - (3). Равновесие, в котором присутствуют только неинфицированные клетки, а инфицированные клетки или частицы вируса всегда существуют, и оно глобально асимптотически устойчиво для R0 ≤ 1 и нестабильно для R0> 1. Нетривиальное равновесие с неинфицированными клетками, инфицированными клетками и вирусными частицами существует для R0 > 1 и локально устойчиво, когда оно существует.
Теперь используем моделирование с реалистичными значениями параметров для дальнейшего изучения вероятного поведения модели. Потребовалось два набора значений параметров для моделирования. Один набор должен был использоваться для получения значения R0 <1, а другой - для R0> 1. Следовательно, значения параметров были получены из разных источников, так что они были бы реалистичными, но их не нужно было существенно изменять по сравнению с опубликованными значениями.
Наши значения параметров были взяты из опубликованной литературы. λ, β, δ, u, a были взяты из [5], а N - из [6]. Значения параметров для рисунка 2: β = 1 × 10–10 сут – 1 дм3, λ = 107 сут – 1 дм – 3, N = 100, γ = 0,2, η = 0,2, a = 0,5 сут – 1, δ = 0,1 сут -1 и u = 5 сут -1