Опыт моделирования динамики профзаболеваний

Г.В.Федорович Д.ф.-м.н., технический директор ООО «НТМ-Защита», Москва

Введение.

В работе предлагается дескриптивная математическая модель развития профессиональных заболеваний (далее ПЗ). Как всякая модель – это упрощенная картина реального явления, используемая для изучения его ключевых свойств. Она представляет собой множество взаимосвязанных предположений об объекте моделирования и отражает некоторые, хотя и не все, свойства реального явления. Цели моделирования:

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

Непосредственным поводом для разработки модели явился вполне результативный опыт моделирования заболеваний с временной утратой трудоспособности (далее ЗВУТ) в работе [1]. В качестве демонстрации возможностей моделирования, отметим, что в логике модели [1] потребовалось пересчитать рутинные данные о ЗВУТ (получаемые в соответствии со специальными методическими рекомендациями углубленного изучения ЗВУТ [2]) во входные параметры модели – вероятности заболеть α и выздороветь β. Эти данные позволили определить выходные параметры – вероятности найти в коллективе здорового Р0 и больного Р1 работника. В ходе построения модели были выявлены важные параметры динамики ЗВУТ - средние длительности цикла «заболевание-выздоровление» L и продолжительности болезни l , с которыми связаны вероятности α и β.

Непосредственная проверка справедливости результатов моделирования на конкретном материале эпидемиологических исследований свидетельствует об адекватности моделирования ЗВУТ при определении профессиональных рисков, связанных с производственной деятельностью [1].
В работе [3] был предложен другой способ описания статистического распределения частоты ПЗ работников по времени контакта с вредными производственными факторами (далее ВПФ). Использовалась модель случайных блуждания вдоль оси Х, по которой откладывался уровень накопления функциональных изменений в организме работников, подвергающихся неблагоприятному воздействию ВПФ. Случайные блуждания, анализируемые методами теории цепей Маркова, приводят к рутинному биометрическому описанию распределения вероятности ПЗ в зависимости от стажа работы с ВПФ. Это открывает возможность определять такие практически важные характеристики, как общий и постажевый коэффициенты профзаболеваемости, ожидаемую продолжительность работы для лиц, доработавших до заданного стажа τ и т.д.
Заметим, что предложенная в [3] модель случайных блужданий была использована , по существу, лишь для аппроксимации статистических данных о распределении частоты ПЗ работников по времени контакта с ВПФ, тем не менее она позволяет оценить качество и достоверность самих исходных данных.

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

1. Динамика ПЗ.

1.1. Целесообразно исходить из вполне очевидного утверждения: профзаболевание – это в первую очередь заболевание и к нему следует относить все, что медицина накопила для описания заболеваний. Например, определения различных стадий болезни.
Общая патология уже давно утвердилась в представлениях о развитии заболевания задолго до его клинического проявления, показав, что доступные клинической диагностике морфологические изменения означают не начало заболевания, а уже ту его стадию, с которой патологические изменения могут превращаться в самоподдерживающийся процесс безотносительно к воздействию этиологического фактора. Началом заболевания следует считать начало воздействия на организм болезнетворного агента. Применительно к проблемам профэпидемиологии можно утверждать, что заболевание начинается при первом контакте работника с ВПФ. Первый период развития ПЗ – субклиническая форма заболевания, протекающая скрытно и не проявляющаяся как физиологическое расстройство. Благодаря защитным приспособительным реакциям организма, у него вначале имеется достаточно ресурсов для подавления возможных проявлений физиологических расстройств.

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

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

В ряде работ, посвященных динамике ПЗ, указывается на последовательные стадии развития заболеваний. Например, в работе [4] представлена характеристика стадий дизадаптоза на примере горнорабочих виброопасных профессий. В [5] отмечается, что результаты выполненных исследований подтверждают стадийностьформирования адаптивно-защитных реакций организма при продолжающемся воздействии факторов рабочей среды. Наличие продромальной стадии ПЗ выявлено у работников нефтехимических производств [6]. Метаболические изменения на клеточном и субклеточном уровне выявлялись при лабораторных исследованиях у лиц без клинических проявлений заболеваний, что свидетельствует о развитии патологических процессов в организме работников.

Из материалов исследования ЗВУТ работников нефтеперерабатывающей промышленности [7] и [8] можно извлечь информацию о количественных характеристиках динамики развития ЗВУТ.

1.2. Как отмечалось во Введении, важный параметр динамики ЗВУТ – средняя длительность цикла «заболевание-выздоровление» L может быть определена по результатам углубленного изучения ЗВУТ – количеству К заболеваний в коллективе из N работников в течение определенного срока Y : L = YN/K . Обратная величина v = 1/L = K/NY имеет смысл частоты ЗВУТ. В работе [7] было обнаружено, что длительность циклов L сокращается (соответственно, частота v растет) у работников по мере увеличения стажа работы во вредных условиях. Впоследствии это было подтверждено в работе [8]. Закон изменения L и v не зависит от модели динамики ПЗ, он должен выполняться во всех моделях. Определим его количественные характеристики.

Для анализа процесса роста частоты ЗВУТ введем счетчик j циклов, произошедших от начала работы до стажа τ. Приведенные в [7] и [8] данные о ЗВУТ относятся к начальному стажу (τ 0 ≤ 3 года) и к стажу τ 1 = 12 лет. Для нефтяников Сибири вначале K0 = 41,8 случаев ЗВУТ (результаты приводятся к численности коллектива N = 100 чел. и длительности наблюдений Y = 1 год). Это соответствует L 0 = 2,4 года (v0 = 0,42 год-1). Через τ1 = 12 лет наблюдалось K1 = 168,3 случая ЗВУТ . Эти 12 лет были заполнены некоторым количеством J циклов, К концу срока период отдельного цикла стал равным L J = 0,59 лет (vJ = 1,7 год-1).
Зависимость v от j может быть любой, лишь бы v росла с ростом j . Однако, имея в распоряжении только две точки функции vj , можно определить параметры только простейшей (линейной) зависимости vj = v0*(1+j*χ) . Примем по-необходимости, эту гипотезу и определим коэффициент χ и количество циклов i , необходимых для заполнения стажа τ 1 = 12 лет. Имеем систему уравнений:

L 0 + L 0 /(1+χ) + L 0 /(1+2*χ) + … + L 0 /(1+J*χ) = τ1 (1)

L 0 /(1+J*χ) = L J (2)

2. Модель.

2.1. При моделировании ЗВУТ в [1] считалось, что есть некоторая переменная Х, характеризующая физиологическое расстройство какой-либо системы жизнедеятельности, которая может принимать два значения: Х = 0, соответствующее отсутствию расстройств и Х = 1, соответствующее наличию расстройства (заболеванию).
При моделировании развития ПЗ в работе в работе [3] переменная Х могла меняться вдоль неограниченной оси со счетным числом состояний. Допускались переходы в соседние состояния.
Объединим эти модели.
Будем считать, что смена уровня расстройства физиологического состояния происходит в момент очередного заболевания. При последующем выздоровлении нормализуются внешние признаки, однако физиологическое состояние организма не возвращается на прежний уровень – появляются патологические изменения, хотя болезнь не проявляется в показателях, характеризующих функциональные изменения в организме. Следующий цикл «заболевание-выздоровление» фиксирует новый уровень патологических изменений в организме - расстройства физиологического состояния. И так далее. Формально этот процесс нарастания физиологических расстройств можно характеризовать введенным выше счетчиком j , который в каждом цикле «заболевание-выздоровление» прирастает на единицу. Таким образом, уровень j – это просто количество циклов, пережитых организмом до настоящего времени. Его необходимо связать с внешними проявлениями скрытых (до поры до времени) физиологических расстройств в организме. В качестве такого параметра, доступного для внешнего наблюдения, предлагается длительность L цикла «заболевание-выздоровление». Действительно, (см. обсуждение выше), по мере увеличения стажа работы длительность цикла сокращается. Это можно использовать для характеристики действия адаптивных механизмов организма, позволяющих ему выдерживать вредные внешние воздействия (вначале успешнее, затем хуже) без видимых последствий.

2.2. Перейдем к математической формулировке модели. Обозначим через pj вероятность найти организм на j-том уровне. Будем считать, что вероятность перехода на уровень j+1 за малый интервал времени (τ,τ+dτ) равна vj*dτ, а вероятность сохранения прежнего уровня равна 1 - vj*dτ. В начальный момент времени (τ=0) организм находится на нулевом уровне, так что p0= 1 , а для j >0 pj = 0 .
Вероятность найти организм на нулевом уровне убывает с течением времени за счет перехода на уровень j = 1. Этот процесс описывается уравнением

dp0/dτ = - v0*p0 (3)

Для уровней с j ≥ 1 изменения вероятностей при переходах с уровня на уровень описываются уравнениями

dpj/dτ = - vj*pj + vj-1*pj-1 (4)

В уравнениях (3) и (4) величины 1/vj определяют характерное время пребывания организма на j-том уровне. Учитывая определение номера уровня j (счетчик количества циклов), можно утверждать, что величина 1/vj равна длительности Lj очередного цикла «заболевание-выздоровление». В отличие от работы [1] , в которой рассматривались циклы постоянной длительности, в предлагаемой модели будем учитывать возможные изменения длительности циклов при переходе с одного уровня физиологических расстройств на следующий (см.выше п.1.2). В рамках простейшего предположения о линейной зависимости vj от j имеем:

vj = v*(1 + χ*j) или L j = L /(1 + χ*j) (5)

Здесь под L = 1/v подразумевается длительность цикла «заболевание-выздоровление» до начала ПЗ (до первых контактов работника с ВПФ). Коэффициент χ в соотношениях (5) это и есть та физиологическая характеристика уровня воздействия ВПФ, с помощью которой осуществляется переход от формального номера состояния j к реальному параметру L j , отражающему адаптивные способности организма. Введение в модель этого коэффициента дает возможность объединить в ней основные функциональные комплексы (внешнюю среду и организм).
2.3. Рассмотрим подробнее уравнения модели с линейной зависимостью vj от j . Для j = 0 уравнение остается тем же (3) с заменой v0 на v, а для j > 0 уравнения (4) имеют вид:

dpj/dτ = - v*(1 + j*χ)*pj + v*[ 1 + (j-1)*χ ]*pj-1 (6)

Решение уравнений (3) и (6) с указанными выше начальными условиями имеет вид

pj(τ) = λj*exp(-vτ )*[1 - exp(-vχτ )] j (7)

где, очевидно, λj0 = 1, а для j ≥ 1 коэффициенты λj удовлетворяют рекуррентным соотношениям

χ*j*λj = [ 1 +χ(j – 1) ]*λj-1 (8)

Распределения вероятностей (7) описывают изменения вероятности пребывания отдельного работника на различных уровнях j физиологического состояния по мере увеличения стажа работы τ . Для статистического описания распределений в коллективе работников по стажу τ , необходимо просуммировать распределения (7) по индексам физиологического состояния j . Количество возможных уровней должно быть ограничено сверху некоторым индексом j0 , соответствующим такому состоянию физиологических расстройств в организме, при которых дальнейшая работа становится невозможной. Работник признается инвалидом и либо уходит на пенсию, либо переводится на другую работу, на которой он меньше подвержен действию ВПФ.
В результате суммирования слагаемых (7) в пределах 0 ≤j ≤ j0 получим функцию распределения n0(τ) работников по стажу, которую, по аналогии с биометрической функцией дожития в демографии (см. напр. [9]), можно назвать функцией дорабатывания.
Непосредственный подсчет суммы слагаемых (7) по индексам j затруднителен из-за сложной зависимости от j коэффициентов λj . Можно, однако, воспользоваться тем, что функция n0(τ) является решением достаточно простого дифференциального уравнения. Для его вывода просуммируем обе стороны уравнений (6) для всех j и прибавим к результату уравнение (3). Получим:

dn0/dτ = - υ*( 1 + j*χ )*pj(τ) (9)

Здесь и ниже следует принять j = j0 . Так как решения pj(τ) для всех j известны (см. (7)), решение n0(τ) легко получается простым интегрированием (5), при соблюдении естественного требования n0(τ→∞) → 0 . Получим:

n0 (τ)=υ*(1+j* )*λjτexp(-υt)*[1-exp(-υχt)]j dt (10)

Интеграл в этом соотношении можно выразить через специальную функцию - неполную бета-функцию Bх(a,b) (см. напр. [10]).Так как при τ= 0 все pj , кроме р0 , равны нулю, а р0 = 1, то должно быть n0(τ=0) =1 . Это обстоятельство дает возможность представитьλj через В1(1/χ, j+1):

(11)

Используя рекуррентные соотношения для бета-функции, можно непосредственно убедиться в выполнении условий (8) для λj. Окончательно имеем

где обозначено x = exp(-υχτ).

Отметим, что неполная бета-функция часто встречается в статистике. Например, интегральные функции вероятности биномиального распределения, F-распределения и распределения Стьюдента выражаются через неполную бета-функцию.
Если функцию распределения n0(τ) рассматривать как аналог биометрической функции дожития s(t) в демографии, то ее производную по времени (9), так же по аналогии с биометрической функцией скорости гибели w(t) =- ds/dt , можно считать скоростью ухода с работы в связи с ПЗ.

Рис.1 Функции дорабатывания nj(t) и скорость ухода с работы dnj/dτ (для j = j0).

Обе функции представлены на графике рис.1. Видно, что статистическое описание динамики развития ПЗ внешне аналогично описанию процессов дожития и смертности в демографии. Для модельного описания ПЗ достаточно нескольких параметров (L , χ , j0), причем их смысл и направления изменения вполне очевидны. Например, величины L и j0 наблюдаемы непосредственно, ухудшение условий труда отражается в модели ростом коэффициента χ и т.д.

3. Свойства модели.

Несмотря на то, что предложенная модель имеет описательный характер, она предоставляет удобные инструменты для проведения практически значимых расчетов. Например, используемые для описания динамики ЗВУТ параметры - средняя длительность цикла «заболевание-выздоровление» L и продолжительность болезни l , непосредственно определяют нетто-тариф взносов δ, страхующих оплату больничных листов заболевшему работнику.
Исходные данные для таких расчетов получаются в результате углубленного изучения заболеваемости с временной утратой трудоспособности, проводимого по методике [2]. Методика определяет тип т.н. «продольных исследований». Регистрируются количество К случаев ЗВУТ и их общая длительность D в коллективе (нормируется на N=100 работников) в течение определенного срока Y (обычно – одного года). Для расчета страховых взносов необходимы результаты типа «поперечных исследований»: длительность заболевания l =D/K и период цикла «заболевание-выздоровление» L = YN/K . Последние данные определяют вероятностные характеристики ЗВУТ: вероятность заболеть в единицу времени a=1/(L+l) , вероятность выздороветь в единицу времени b = 1/l и , предполагая Марковский характер процесса, вероятность найти здорового P0 = β/(α+β) или больного P1 = α/(α+β) работника в коллективе. Если через q обозначить заработок (начисляемый) отдельного работающего, то суммарный заработок коллектива из N работников будет q*N*P0 . После вычитания страховой доли δ зарплата всего коллектива будет (1-δ)*q*N . Из условия баланса заработка и зарплаты получим : δ = 1 – Р0 = Р1 .
В середине 90-х годов прошлого века специалистами НИИ Медицины труда РАМН были проведены массовые и длительные по срокам наблюдения статистические исследования роли санитарно-гигиенических условий труда в формировании уровней заболеваемости с временной утратой трудоспособности [11]. На основе уровней вредных производственных факторов непосредственно на рабочих местах, определяемых с помощью инструментальных замеров, устанавливался класс вредности условий труда и, в сопоставлении с показателями здоровья и величины утраты трудоспособности, определялась величина вероятности для работника негативных последствий для здоровья.
Диады {K,D} были определены для различных производств с различными классами условий труда (КУТ). Эти данные из [11] сведены в таблицу 1. В последней строке – рассчитанный по описанной выше методике нетто-тариф взносов.

Таблица 1

Рассчитанные отчисления в среднем близки к реальным взносам 2,9 % в ФСС. Результаты, приведенные в табл.1, могут служить исходными данными для оценки финансовой выгоды (или потерь) от мер по улучшению условий труда на производстве и возможного времени их окупаемости.

2.2.Как показано выше (см. п.1.4), рассматриваемая здесь модель приводит к стандартному биометрическому описанию распределения работников по стажу. В предыдущей статье [3] было показано, как, используя такое описание (фактически - функцию дорабатывания n0(τ)) , можно проводить оценку страховых нетто-тарифов, покрывающих пенсионные выплаты работникам, вышедшим на пенсию из-за хронических ПЗ. Такие оценки можно воспроизвести, используя соотношения (10) и (12) для функции n0(τ) . Необходимо, однако, знать индекс j0 , определяющий уровень физиологических расстройств в организме, при которых идентифицируется ПЗ и дальнейшая работа становится невозможной. Для этого нужна информация, внешняя по отношению к модели. Например, можно использовать данные из сборника [12] о распределении ПЗ по времени контакта (стажу) работников с ВПФ. Они приведены к однолетнему интервалу времени, т.е. отображают именно скорость (9) ухода с работы в связи с ПЗ. В графическом виде результаты представлены в работе [1]. Анализ кривых, интерполирующих скорость, показывает, что при работе во вредных условиях труда (КУТ 3.1-4) максимум приходится на время контакта (стаж) ≈ 25 лет. При модельном описании динамики ПЗ этому стажу соответствует положение максимума функции pj(τ) (см. рис.1). Непосредственное дифференцирование функции (7) по времени τ определяет положение τmax максимума для вероятностей pj(τ) на различных уровнях j.

1 + j* χ = exp( v*χ* τmax ) (13)

Используя это соотношение, можно определить одну из входящих в него величин, если известны остальные. Например, для полученных выше значений v = 0.3 лет-1, χ = 0.33 , τm ≈ 25 , из (13) следует j0 ≈ 34 . Таким образом, можно ожидать, что после того, как работник пережил 30 – 40 случаев ЗВУТ, у него с большой вероятностью будет установлено ПЗ. Длительность циклов «заболевание-выздоровление» при этом сократится до [v*(1+ χ*j0)]-1 ≈ 0.27 года = 100 дней, т.е. работник будет болеть не меньше 3-х раз в год. Это может иметь важное прогностическое значение как сигнализатор опасности профессионально обусловленной инвалидности на ранних стадиях развития ПЗ.

Заключение.

В качестве основного результата проведенного рассмотрения возможностей моделирования динамики развития ПЗ, можно считать демонстрацию плодотворности взгляда на ПЗ как на единый процесс – от первого контакта работника с ВПФ до выхода на пенсию или инвалидности.
При таком подходе выявляется важнейший параметр, характеризующий динамику развития ПЗ – длительность циклов «заболевание-выздоровление». Их укорочение свидетельствует о развитии ПЗ, несмотря на возможное отсутствие (вначале) внешних проявлений серьезных функциональных нарушений в организме.
Выявление вариаций длительности циклов «заболевание-выздоровление» в зависимости от их порядкового номера j дает важную прогностическую характеристику динамики ПЗ, сигнализирующую об опасности профессионально обусловленной инвалидности на ранних стадиях ее развития. Количественная характеристика этой зависимости – коэффициент χ линейной связи частоты ЗВУТ с номером j указывает на степень влияния ВПФ.
Несомненно необходимо нозологическое структурирование исходных статистических данные, в том числе и полученных методом углубленного изучения ЗВУТ.

Литература.

  1. Федорович Г.В. Аппарат теории цепей Маркова в профэпидемиологии // БиОТ. – 2013. – № 1. – С. 61 –65.
  2. Догле Н.Ф., Юркевич А.Я. Углубленное изучение заболеваемости с временной утратой трудоспособности и выявление роли санитарно-гигиенических условий труда в формировании ее уровней.- М.: Медицина, 1984 - 256 с.
  3. Федорович Г.В. Основания актуарных расчетов рисков профессиональных заболеваний // БиОТ. – 2013. – № 2. (в печати).
  4. Коневских Л.А., Оранский И.Е., Макогон И.С. Адаптивные возможности сердечно-сосудистой системы у горнорабочих виброопасных профессий // Медицина труда и промышленная экология – 2013. – №2. – С. 32 – 37.
  5. Юдина Т.В., Сааркоппель Л.М., Крючкова Е.Н. Иммунореактивность организма рабочих при производстве цемента // Медицина труда и промышленная экология – 2013. – №3. – С. 6 – 11.
  6. Тимашева Г.В., Кузьмина Л.П., Каримова Л.К., Бадамшина Г.Г. Роль лабораторных исследований в диагностике ранних метаболических нарушений у работников нефтехимического производства // Медицина труда и промышленная экология – 2013. – №3. – С. 15 – 20.
  7. Овчаров Е.А. Характеристика заболеваемости с временной утратой трудоспособности нефтяников Западной Сибири //Здравоохр. Рос. Федерации. - 1996. - № 5. - С. 35-38.
  8. Даткаева Г.М., Булешов М. А.,Шаимерданова Б.Е. Динамика и структура заболеваемости с временной утратой трудоспособности работников нефтеперерабатывающей промышленности Южного Казахстана //
  9. Сакович В.А., Смирнова О.А. Математическое моделирование влияния радиации на продолжительность жизни млекопитающих // Физика элементарных частиц и атомного ядра. – 2003. – Т.34. – Вып.6. – С.1436 -1481.
  10. Градштейн И.С., Рыжик И.М. Таблицы интегралов, сумм, рядов и произведений.-М.: Физматгиз, 1963 – 1100 с.
  11. Молодкина Н.Н., Радионова Г.И., Денисов Э.И. Обоснование критериев профессионального риска. В кн. Измеров Н.Ф. (ред) Профессиональный риск. – М.: Социздат, 2001. – С. 48 - 55.
  12. Верещагин А.И.(ред).Профессиональные заболевания и их распределение по классам условий труда в Российской Федерации в 2009 году / Информационный сборник статистических материалов. - М.: Роспотребнадзор. - 2010. -108 с.