вторник, 14 апреля 2020 г.

Модель SIR и эпидемия COVID-19: продолжение

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

График ниже будет обновляться в соответствии с новыми данными, так что имеет смысл поглядывать время от времени и в дальнейшем. Текущее обновление, с уточненной сегментацией этапов эпидемии и ее математических моделей: 27 апреля. 

По состоянию на 2 мая к описанию сделано добавление, представляющее анализ ситуации в связи с новыми данными, добавившимися за время, прошедшее с первого варианта описания модели.


Итак, начнем снова с SIR-модели в представлении уравнения Ферхюльста, которая для кумулятивного числа инфицированных имеет вид
$$\frac{d\Sigma I}{dt}=r\Sigma I\left(1-\frac{\mathbf{R}}{K}\right)$$
 с  решением 
$$\Sigma I(t)=\frac{K}{1+e^{-r(t-t_m)}}.$$

Однако, для поставленной цели это решение все же удобнее переписать в виде, выраженном через начальное число инфицированных в качестве одной из констант:
$$\Sigma I(t)=\frac{\Sigma I(0)e^{rt}}{1+\frac{\Sigma I(0)}{K}\left(e^{rt}-1\right)}$$
так как его логарифм имеет вид
$$\ln(\Sigma I(t))=
\ln(\Sigma I(0))+rt-\ln\left(1+\frac{d\Sigma I(0)}{K}\left(e^{rt}-1\right)  \right).$$

Первое слагаемое - постоянное, второе - линейная функция от времени (что должно быть хорошо заметно на графике в полулогарифмических координатах), третье, которое вычитается, приводит к уменьшению $\ln(d\Sigma I(t))$ при приближении решения к насыщению. При этом когда $rt<<1$ или же, более обще, $\frac{d\Sigma I(0)}{K}\left(e^{rt}-1\right) <<1$, это третье слагаемое близко к нулю так как $\ln(1)=0$, т.е. прямая на полулогаримфическом графике как раз свидетельствует о том, что динамика роста все находится на начальном этапе - этапе экспоненциального роста, относительно вдали от максимума текущего числа инфицированных и, тем более, от конца эпидемического процесса, соответствующего насыщению
$\Sigma I=\mathrm{max}{(\Sigma I)}$.

Теперь непосредственно к данным, касающимся, России - для вычислений, соответствующих уравнениям выше, и построения графика, скрипт для MATLAB из предыдущего материала заменяется следующим (изменения связаны с обновлением данных в ECDC 20 апреля 2020 - были заполнены пропуски строк в сплошной серии дней - в дни, в которые данные от национальных служб здравоохранения данных не поступало, указаны нули вместо пропущеных строк) - он считывает даты из файла ECDC и строит графики, а также проводит их линейный фит и находит его наклон для выделенных временных участков (в новом варианте искомые границы временных интервалов задаются напрямую указанием дат):

clear;close all;clc;
% Variant valid since 20/04/12 - the gaps in data were eliminated

TBL=readtable('COVID-19-geographic-disbtribution-worldwide.csv');
countries=TBL.geoId;
geoIdcode='RU';
idx=find(strcmp(countries,geoIdcode));
len=length(idx);
I=TBL.cases(idx);
I=flipud(I);
Dates=TBL.dateRep(idx);
Dates=flipud(Dates);
t=0:len-1;
t=t';
D=TBL.day(idx);
D=flipud(D);
U=cumtrapz(t,I);

dat1=find(strcmp(Dates,'27/02/2020'));
dat2=find(strcmp(Dates,'07/03/2020'));
dat2a=find(strcmp(Dates,'15/03/2020'));
dat2b=find(strcmp(Dates,'25/03/2020'));
dat3=find(strcmp(Dates,'04/04/2020'));
dat4=find(strcmp(Dates,'13/04/2020'));
dat5=find(strcmp(Dates,'18/04/2020'));
dat6=find(strcmp(Dates,'24/04/2020'));

    figure
    plot(t(dat1:dat2a),log(U(dat1:dat2a)),'o-')
    hold on
    plot(t(dat2a:dat3),log(U(dat2a:dat3)),'o','color','green')
    p1=polyfit(t(dat2a:dat3),log(U(dat2a:dat3)),1);
    plot(t(dat2:dat4),polyval(p1,t(dat2:dat4)),'-','color','green','LineWidth',1)
    plot(t(dat3:dat4),log(U(dat3:dat4)),'o','color','red')
    p2=polyfit(t(dat3:dat4),log(U(dat3:dat4)),1);
    plot([t(dat1:end);t(end)+3],polyval(p2,[t(dat1:end);t(end)+3]),'-','color','red','LineWidth',1)
    plot(t(dat4:end),log(U(dat4:end)),'o','color','black')
 
    plot([t(dat2b) t(dat2b)],[0 14],':','color','black')

    set(gca,'XTick',[t(dat1), t(dat2a), t(dat3), t(dat5)],...
        'XTickLabel',{Dates{dat1},Dates{dat2a},Dates{dat3},Dates{dat5}})
    xlabel('Date')
    ylabel('\Sigma I')
    title(['Current date: ',cell2mat(TBL.dateRep(idx(1)))])
    set(gca,'FontSize',11)
    grid on

Результат расчета, выполненого 17 апреля (еще используя данные ECDC c пропусками) - на графике ниже.


Явно видна сегментация на три (пока) участка с различным качественным поведением:

  • Синие маркеры, соединенные линией - ускоренный (быстрее чем экспоненциальный) рост на начальном этапе - это нормальный эффект усреднения по массе независимых неупорядоченных данных, в каждом из источников кривая может расти экспоненциально, но из-за стохастического разброса кривых роста их суммирование приводит к квадратичному росту в полулогарифмических координатах (линейна часть зависит от средней скорости роста, квадратичная - от дисперсии скоростей). Соответствующее обоснование в подробностях может быть найдены в статье (она относится к росту другого патогена - микобактерии, вызывающей туберкулез и в лабораторной культуре, но математика популяционной динамики - соверешенно аналогична):
  • Зеленые маркеры от 7 марта до 4 апреля - уже указанный выше стандартный экспоненциальный рост, линейный фит, проведенный через них, дает коэффициент роста $r_1$, указанный на графике. Такая установившаяся зависимость может проистекать из двух причин: а) один из компонентов (Москва?) подавляет все остальные слагаемые; б) скорости роста примерно выровнялись для разных слагаемых, или же их комбинацией. 
  • Красные маркеры, начиная с 5 апреля -  по первым из них можно было решить, что, раз они пошли под зеленую линию, то это как раз то самое замедление процесса, идущего однородно во времени, учитываемое третьим слагаемым в уравнении выше. Но учет дальнейших данных показывает, что красные ход красных маркеров не следует кривой замедленного роста, а ложится опять таки на прямую - см. линейный фит за период 5-13 апреля. Т.е. рост - экспоненциальный с высокой точностью, однако, с заметно меньшей скоростью прироста $r_2<r_1$. Учитывая даты, судя по всему это как раз и есть проявившийся эффект карантинных мер. Почему же в ежедневных сводках виден заметно больших прирост день ото дня по сравнению с периодом до 5 апреля? Это происходит потому, что к 5 апреля уже накопилось достаточно большое число носителей, а рост все равно экспоненциальный (т.е. на простом примере - если (числа просто для примера, абстрактные), при скорости ежедневного роста - коэффициенте контактного заражения - 5 имелось 5 носителей, то они дадут в следующий момент времени $5\times 5=25$ новых, но 25 при меньшем коэффициенте заражения, допустим, 2, дадут за тот же промежуток времени уже $25\times 2=50$). 
  • Черные маркеры -  показывают, следует ли продолжения после 13 апреля линейному фиту, отвечающему экспоненциальному росту, в дальнейшем. Замедления на этом графике будут проявляться, когда эти маркеры пойдут стабильно отклоняться от красной линии вниз. 
К 27 апреля ситуация изменилась - текущий вид этого графика:

Участки с качественно различным поведением, описанные выше обозначены теми же цветами (кстати, видно, что накопилась достаточная статистика, чтобы данные аккуратно легли на быстрый экспоненциальный рост, к 25 марта - зеленые маркеры стали очень мало отклоняться от фитующей прямой). Выявленный слом кривой в полулогарифмическом масштабе, начавшийся 4 апреля и профитованный как красная прямая на интервале от 4 до 13 апреля (красные маркеры) продолжался до примерно 18 апреля, когда стал явно заметным начинающийся загиб вниз, существенно ускоряющийся с тех пор, свидетельствующий о прекращении медленной экспоненциальной фазы.

Поэтому имеет смысл перестроить график уже в линейном масштабе:

Цвета обозначений сохранены и в таком масштабе еще более нагляден сход с экспоненты (красная линия) и переход к линейному росту кумулятивного числа инфицированных (то есть, флуктуирующему вокруг среднего, но более-менее постоянного числа ежедневно выявляемых инфицированных). Соответствующая линейная зависимость профитована прямой на интервале от 18 до 24 апреля (синие маркеры), за последующие три дня два последних маркера чуть выше, но, в силу еще малого количества, в принципе пока не показывают заметных отклонений за пределами нормальных флуктуаций как и у всякой статистической величины.

Т.к. производная постоянной $I=\mathrm{const}$ равна нулю  $dI/dt=0$ (или же в терминах кумулятивного числа $d^2\Sigma I/dt^2=0$, т.е. $\Sigma I=kt+b$), то это свидетельствует об окрестности максимума - но эта окрестность достаточно размазана - с неопределенностью не меньше недели по времени (минимум - если середина линейного интервала, а т.к. данные еще за линейный режим не вышли - коридор расширяется до удвоенного срока наблюдения линейности). Однако, упрощенная модель SIR-Ферхюльста, рассмотренная ранее такую ситуацию не воспроизводит - т.к. максимум достаточно точечный, но при этом такое решение - следствие допущения в инварианте
$$\log(S)+k\tau R=\log(S_0)$$
равенства $S_0\approx 1$ (при нормировке общей восприимчивой популяции на единицу), что соответствует самому началу эпидемии (еще практически нет инфицированных $I<<1$ и переболевших $R=0$) и, при этом параметры распространения постоянны. К моменту изменения скорости экспоненциального роста это, очевидно, не выполняется. Поэтому модель имеет смысл рассматривать стартуя с новых начальных условий, $R'=R-R_{13}$ (это обеспечит $R'=0$ в новый начальный момент) и $S'_0=S_{13}$, существенно не равное единице (число восприимчивых в новый начальные момент $S_{13}=1-I_{13}-R_{13}$ при той же нормировке на единицу).

В этом случае, в точном уравнении SIR-модели
$$
\frac{dR}{dt}==\tau^{-1}\left[1-R'-S'_0e^{-k\tau R'}\right].
$$
при разложении экспоненты в ряд Тейлора единица в скобках не сокращается и интеграл метода разделенных переменных принимает вид
$$
\int\limits_0^R\frac{dR}{(1-S_0)+(S_0\mathbb{R}_0-1)R-S_0\frac{\mathbb{R}_0^2}{2}R^2}=\tau t,
$$
из которой можно выделить полный квадрат,
$$
\int\limits_0^R\frac{dR}{\left[\sqrt{\frac{2(1-S_0)}{S_0\mathbb{R}_0^2}+\left(\frac{S_0\mathbb{R}_0-1}{S_0\mathbb{R}_0^2}\right)^2}\right]^2-\left[R-\frac{S_0\mathbb{R}_0-1}{2S_0\mathbb{R}_0^2}\right]^2}=
\frac{S_0\mathbb{R}_0^2\tau}{2} t,
$$
то есть получить табличный интеграл, равный гиперболическому арктангенсу (обозначив $\xi=R-b$):
$$
\int\frac{d\xi}{a^2+\xi^2}=\frac{1}{a}\tanh^{-1}\left(\frac{\xi}{a}\right).
$$

Соответственно, искомое решение выражается через гиперболический тангенс:
$$
R=b+a\tanh\left[\frac{aS_0\mathbb{R}_0^2\tau}{2} t-\tanh^{-1}\left(\frac{b}{a}\right)\right].
$$

Переходя от $R$ а пропорциональному ему кумулятивному числу инфицированных к текущему моменту,
$$
\Sigma I=\Sigma I_{13}+a\tanh\left[\frac{aS'_0\mathbb{R}_0^2\tau}{2}(t-t'_m)\right],
$$
где $t'm$ - истинный временной момент пика заболевания. При этом, т.к. гиперболический тангенс, аналогично обычному тангенсу и обычному синусу, в достаточно протяженной окрестности нуля линеен, получается воспроизведение в модели искомого свойства линейного роста.

Анализируя данные - т.е. приравнивая $\frac{aS'_0\mathbb{R}_0^2\tau}{2}$ найденному наклону синей фитующей прямой на графике выше, а также подбирая $t'_m$ и $a$ таким образом, чтобы решение через гиперболический тангенс при $t<<t'_m$ переходило в экспоненциальный рост, найденный как красная фитующая экспонента,
получаем решение, показанное на графике выше малиновой штриховой линией, пока что успешно воспроизводящей наблюдаемый процесс.

 Дальнейшие наблюдения покажут, будет ли следовать динамика такому решению.

Добавления от 2 мая.

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



Автор: на