Введение
Модели «доза-эффект» на основе логит- или пробит-функций связи имеют более чем полувековую историю (Finney, 1947, 1971; Беленький, 1963). Развитие этой методологии от использования полуэвристических закономерностей до полноценных регрессионных моделей на основе принципа максимального правдоподобия не столько увеличило точность оценки токсикометрических показателей, сколько создало предпосылки для использования мощных и эффективных способов контроля за воздействием случайных факторов.
При решении задач экотоксикологии невозможно провести «стерильный» эксперимент с химически стандартным токсическим веществом. Трудно также оценить уровень опасного воздействия в условиях видовой, возрастной или гендерной специфичности популяций, неоднородности природной среды, влияния сезонно-климатических изменений и прочих сопутствующих факторов. Современный подход к решению этих проблем связан с использованием регрессионных моделей со смешанными параметрами (mixed effects model – Pinheiro, Bates, 2000; Demidenko, 2013; Четвериков, 2015), которые создают предпосылки к существенному повышению уровня обобщения целевых результатов исследования и их распространению на другие пространственные, временные или биологические уровни организации. Интерес к смешанным моделям постоянно растет: по данным Google Scholar, только в 2017 г. было опубликовано 68 000 статей, в том числе по экологии – 17 800 статей (из них на русском языке – не более 3).
Сущность данного метода заключается в том, что эффекты (факторы), оказывающие влияние на зависимую переменную, условно разделяются на два типа: фиксированные и случайные. Здесь фиксированные эффекты – это то, что обычно является предметом интереса исследователя, то есть те независимые переменные, уровни которых он устанавливает или контролирует (в нашем случае – доза яда гадюки). Для случайных эффектов обычно затруднительно заранее составить в исследовании план варьирования уровней факторов, и они могут принимать случайные значения из числа возможных в генеральной совокупности. Философско-методологические аспекты разделения эффектов на типы весьма неоднозначны (Gelman, 2005), поэтому ниже мы будем интерпретировать как случайные все сопутствующие факторы, в той или иной степени влияющие на процесс интоксикации и гибели подопытных животных (их пол, подвид и район отлова гадюк, цвет яда и т. д.).
Цель работы – оценить результаты применения моделей со смешанными параметрами на конкретных выборках экотоксикологических данных. Основные задачи исследования: во-первых, показать стратегию анализа отдельных вариантов моделей и возможности их интерпретации; во-вторых, проверить надежность и полезность метода для решения отдельных вопросов токсинологии. Содержательные экологические выводы, относящиеся к пространственному распределению отдельных популяций гадюк и их систематическому статусу, мы выносим за рамки настоящей публикации.
Материалы
Образцы ядовитого секрета были получены от обыкновенных гадюк Vipera berus (Linnaeus, 1758), отловленных в 12 разных географических регионах (далее – переменная местооб) европейской части РФ. Подвидовую принадлежность (переменная подвид) гадюк изучаемых популяций – V. b. berus (Linnaeus, 1758) и V. b. nikolskii Vedmederja et al., 1986 – устанавливали по морфологическим признакам и цвету яда (Bakiev et al., 2005; Milto, Zinenko, 2005). Разделение на два подвида следует воспринимать как условное, поскольку часть популяций находится в зоне интерградации этих подвидов и обладает смешанными признаками. Так, в Пензенском районе Пензенской области выявлена популяция, в которой одна часть гадюк имеет желтый яд (образец 8), характерный для номинативной формы V. b. berus, другая – бесцветный яд (образец 7), характерный для лесостепной формы V. b. nikolskii.
Препараты ядовитого секрета вводили внутрибрюшинно белым лабораторным мышам, самцам и самкам (переменная пол) массой 20.0 ± 1.0 г. Для каждого образца яда гадюк, отловленных в разных регионах, было проверено 5 эффективных доз по 5 мышей в каждой группе, в том числе 86 групп для яда V. b. berus в диапазоне доз от 0.5 до 2.5 мкг/г и 40 групп для яда V. b. nikolskii – от 0.3 до 1.5 мкг/г.
Методы
Предварительно силу воздействия яда оценивали по показателю среднесмертельной токсичности LD50, рассчитанному с использованием статистических моделей «доза-эффект» в трех вариантах:
Вариант А использует линейную зависимость, установленную в середине ХХ в.:
Y = b0 + b1 X (Беленький, 1963),
где Y – так называемые рабочие пробиты Блисса (Bliss, 1935), X – натуральное или прологарифмированное значение воздействующей дозы.
Значение LD50 и его ошибку рассчитывали как
LD50 = (5 – b0)/b1 ; SDLD50 = (LD86 – LD16)/(2N)0.5 (Бабич и др., 2003).
Вариант B основан на стандартной обобщенной модели регрессии с пробит-функцией связи (Мастицкий, Шитиков, 2015):
P = Φ(b0 + b1 X) или Φ-1(P) = b0 + b1 X ,
где Р – вероятность гибели животных, Φ() – интегральная функция плотности стандартного нормального распределения N(μ, σ); b0 = - μ / σ; b1 = 1 / σ; Φ-1(P) – обратная функция от Ф, или пробит. В отличие от варианта А, коэффициенты b0 и b1 оценивали методом нахождения максимума правдоподобия относительно непосредственно целевого отклика (т. е. вероятности эффекта Р) в предположении о биномиальном распределении данных.
Ошибку изоэффективной дозы LDP для произвольно заданного значения Р (т. е. Р = 0.5 при нахождении LD50) оценивали дельта-методом:
SDLD50 = (Jт V J)0.5,
где J – вектор производных функции регрессии относительно вектора параметров b; для рассматриваемой модели J = [1/b1 ; X/b1]; V – дисперсионно-ковариационная матрица коэффициентов логистической модели.
Доверительные интервалы LD50 рассчитывали по обычной формуле
CILD50 = LD50 ± tα/2,f SDLD50,
где tα/2,f – квантиль t-распределения при α/2 и f степенях свободы.
Вариант C также основан на модели логистической регрессии, но доверительные интервалы CILD50 оценивались с использованием процедуры, описанной Д. Финни (Finney, 1971, p. 76–79; Robertson et al., 2003).
В приложении представлены скрипты на языке R для всех трех вариантов функций расчета LD50, оснащенные подробными комментариями. Отметим, что в отличие от функции Probit_Bel(), выполняющей расчеты по Н. Л. Беленькому, на вход остальных двух функций можно подать вектор с любым набором задаваемых вероятностей эффекта, что удобно для построения графиков. Так при p = seq(1,99,1) функции LD_glm() и LD_Rob() могут выполнить расчет таблицы из 99 изоэффективных доз от 1 до 99 %.
Для простой схемы группировки по уровням линейную модель LMEM (linear mixed-effects model) со смешанными эффектами можно представить в матричной форме как:
y = X b + Z b + ε, b ~ N(0, ψ) , ε ~ N(0, Λσ2),
где y – N-мерный вектор зависимой переменной, X – матрица фиксированных переменных; b – вектор параметров модели для фиксированных эффектов; Z – матрица, состоящая из нулей и единиц и описывающая группировку случайных факторов по уровням в соответствии с планом эксперимента, b – коэффициенты модели, относящиеся к случайным эффектам. Оценка значимости параметров модели и ранжирование случайных эффектов по степени их выраженности выполнялись с использованием традиционных показателей качества подгонки обобщенных моделей – информационного критерия Акаике AIC и девиансы (deviance) D, которая является обобщением остаточной суммы квадратов при использовании метода максимального правдоподобия.
Построение моделей со случайными эффектами выполнялось с использованием пакетов lme4 и MuMIn статистической среды R вер. 3.03. Представленный скрипт содержит функции на алгоритмическом языке R для расчета LD50, других изоэффективных доз и их доверительных интервалов тремя различными методами.
Результаты
Расчет значений LD50, представленных в табл. 1 отдельно по V. b. berus и V. b. nikolskii для всех групп местообитаний, показал наличие существенных различий в силе действия ядовитого секрета по подвидам. При этом значения норматива, полученные различными вариантами расчета, практически совпадали.
Таблица 1. Значения LD50, полученные различными способами расчета для двух подвидов гадюк (SE – стандартная ошибка LD50, LCL и UCL – нижняя и верхняя границы 95 % доверительного интервала)
| Подвид | Вариант | LD50 | SE | LCL | UCL |
| Vipera berus berus | А | 1.540 | 0.0832 | 1.38 | 1.70 |
| B | 1.576 | 0.0542 | 1.47 | 1.68 | |
| C | 1.576 | 1.47 | 1.69 | ||
| Vipera berus nikolskii | А | 0.967 | 0.0588 | 0.85 | 1.08 |
| B | 0.972 | 0.0451 | 0.88 | 1.06 | |
| C | 0.972 | 0.88 | 1.07 |
На основании такого точечного показателя, как LD50, нельзя получить достаточно полное представление о сравнительной силе действия ядов, поскольку важно учитывать также скорость развития эффекта при увеличении дозы. На рис. 1 показаны сигмоидальные кривые «доза-эффект», где видно, что на уровне 16 % смертности различия между ядовитыми секретами обоих подвидов минимальны, тогда как 84 % изоэффективные дозы отличаются весьма значительно.
Рис. 1. Кривые «доза-эффект» и их 95 % доверительные интервалы для популяций с преобладанием диагностических признаков V. b. berus или V. b. nikolskii, построенные с использованием логистической регрессии
Fig. 1. "Dose-response" curves and their 95 % confidence intervals for populations with a predominance of diagnostic features of V. b. berus or V. b. nikolskii, built using logistic regression
Далее были рассмотрены вопросы о том, являются ли статистически значимыми полученные зависимости «доза-эффект» и насколько сильно сопутствующие факторы могут влиять на их характер. В табл. 2 приведены результаты сравнительного статистического анализа четырех основных моделей по возрастанию их параметричности. Как принято в R, структура регрессионных моделей задавалась в виде «формулы», в которой зависимая переменная (pd – доля погибших животных) отделяется от правой части знаком ~, а предикторы разделяются знаком +.
Таблица 2. Оценка статистической значимости включения в модель случайных эффектов
| Модель | Формула | AIC | Девианса D | χ2 | P(>χ2) |
| m0 | pd ~ 1 | 587.08 | 437.84 | – | - |
| m1 | pd ~ доза | 312.58 | 161.35 | 276.49 | ≈ 0 |
| m2a | pd ~ доза + подвид | 272.14 | 118.9 | 42.45 | ≈ 0 |
| m2b | pd ~ подвид + подвид:доза – 1 | 264.17 | 108.93 | 9.97 | 0.0016 |
В представленном ряду значения AIC-критерия постоянно уменьшаются, что свидетельствует о поэтапном улучшении модели. Статистическую значимость p включения каждого нового члена (вначале доза, а затем вид) оценивали по разности остаточных девианс (Di-1 – Di). Последняя имеет распределение χ2 и, если она достаточно велика (p < 0.05), то доля вариации, объясняемая новым членом модели, считалась достоверно существенной. На рис. 2 показано, как случайный эффект вид проявляется как статистически значимый параллельный сдвиг линий «доза-эффект» за счет коррекции свободного члена в модели m2a. Модель m2b дополнительно учитывает различия в скорости развития эффекта при увеличении дозы за счет коррекции коэффициента угла наклона. Сравнение моделей m1 и m2a проверяет гипотезу «параллельности», а m1 и m2b – гипотезу «эквивалентности» токсических процессов (Robertson et al., 2003).
Рис. 2. Графики зависимости пробит-значений вероятности смерти мышей при действии яда гадюк V. b. berus или V. b. nikolskii. Пунктиром показаны зависимости, построенные в предположении m2a, что подвидовые различия сводятся к параллельному сдвигу на величину разности LD50. Сплошные линии соответствуют модели m2b, учитывающей дополнительно различия в угле наклона
Fig. 2. Graphs of dependence of probit values of the mice death probability under the action of V. b. berus or V. b. nikolskii venoms. The dotted line shows the dependencies constructed under the m2a supposition that the subspecies differences are reduced to parallel shift by the magnitude of the DL50 difference. The solid lines correspond to the m2b model additionally considering the differences in the slope angle
Смешанную модель оптимальной структуры получили методом «всех возможных регрессий», выполнив полный перебор комбинаций случайных факторов. Модели-претенденты ранжировали по уменьшению скорректированного AICc-критерия (табл. 3). Относительное правдоподобие (strength of evidence (Burnham, Anderson, 2002)) каждой модели mmi по отношению к лучшей модели оценивали как li = exp(-0.5 Δi), где Δi = (AICci – AICcmin). На основании этого можно сделать вывод, что модель mm1, учитывающая эффекты вид и местооб, в 2.2 раза лучше объясняет результаты эксперимента, чем модель mm2, и в 7×106 раза более обоснована (с точки зрения соотношения вероятностей соответствующих гипотез), чем продемонстрированная на рис. 2 модель mm7.
Таблица 3. Ранжирование вариантов моделей со смешанными эффектами по уменьшению информационного критерия Акаике AICc
| Модель | Св. член | доза | Случайные эффекты | Степеней свободы df | Критерий AICc | Δ | |||
| подвид | местооб | пол | |||||||
| mm1 | -2.39 | 1.964 | + | + | 8 | 262.7 | 0 | ||
| mm2 | -2.407 | 1.863 | + | 5 | 263.5 | 0.78 | |||
| mm3 | -2.576 | + | + | 7 | 266.9 | 4.16 | |||
| mm4 | -2.412 | 1.825 | + | + | 8 | 268.6 | 5.85 | ||
| mm5 | -2.391 | 1.959 | + | + | + | 11 | 269.7 | 6.97 | |
| mm6 | -2.429 | + | + | + | 10 | 273.9 | 11.16 | ||
| mm7 | -2.195 | 1.825 | + | 5 | 276.2 | 13.46 | |||
Интерпретация модели со смешанными эффектами основана на том, что для каждого уровня случайного фактора, отмеченного знаком + в табл. 3, подбирается система коэффициентов, осуществляющих коррекцию модели с фиксированными параметрами. Например, фиксированная часть модели mm1 имеет вид Φ-1(P) = -2.39 + 1.964·доза. Если использовали яд V. b. berus, то к этому уравнению добавлялся член (-0.132 – 0.313·доза), а если V. b. nikolskii, то (0.129 + 0.305·доза). Аналогично, для уровня Аткарск фактора местооб – образец 2 – добавляется член (-0.1·доза), для Самара – образец 10 – (+0.32·доза) и т. д. И тогда полное уравнение регрессии с учетом случайных параметров для комбинации «вероятность смерти от яда гадюк V. b. nikolskii, отловленных в Актарске» приобретает вид
Φ-1(P) = (0.129 - 2.39) + (1.964 + 0.305 – 0.1)·доза.
Семейство таких зависимостей «доза-эффект» по результатам проведенного эксперимента представлено на рис. 3.
Рис. 3. Зависимости «доза-эффект» для яда обыкновенных гадюк из популяций, условно отнесенных к подвидам V. b. berus или V. b. nikolskii
Fig. 3. «Dose-response» dependences for the adder venom from populations conventionally related to V. b. berus or V. b. nikolskii subspecies
Пол подопытных животных статистически значимого влияния на развитие токсического процесса не оказывает и в окончательно отобранные модели не включается.
Обсуждение
Анализ закономерностей, полученных с помощью модели mm1, позволяет подтвердить ряд содержательных гипотез по пространственному распределению отдельных популяций гадюк и их систематическому статусу. В табл. 1 показано статистически значимое различие величин LD50 для подвидов гадюк V. b. berus и V. b. nikolskii. Однако на самом деле анализируемые выборки частично сделаны из обширной зоны интерградации этих двух форм, и определить систематический статус отдельных популяций можно было в значительной мере условно. В усилении токсичности ядовитого секрета по отношению к мышам выявляется географический тренд (см. рис. 3), имеющий направление от Пермского края к Липецкой области. Он хорошо отражает современные представления об истории расселения и гибридизации более древней гадюки Никольского с относительно молодой номинативной формой обыкновенной гадюки, отличающейся менее токсичным для мышей ядом.
Таксономическую принадлежность форм со смешанными признаками berus и nikolskii предстоит еще уточнить для многих местообитаний, поскольку сами диагностические признаки последней формы нуждаются в уточнении. В частности, к настоящему времени стали известны популяции V. b. nikolskii без меланистов (Zinenko et al., 2010); значит, черная окраска взрослых особей не может являться диагностическим признаком гадюки Никольского. Кроме того, ее типовые экземпляры добыты в зоне контакта V. b. berus и V. b. nikolskii (Bakiev et al., 2005; Milto, Zinenko, 2005). В связи с возможной перспективой выявления более «чистых» популяций лесостепной формы к западу от типовой территории гадюки Никольского (Zinenko et al., 2010) следует обратить внимание на ст. 23.8 Международного кодекса зоологической номенклатуры (2004), где указывается, что название видовой группы, впоследствии оказавшейся гибридом, не должно употребляться как валидное название ни для одного из родительских видов. Поэтому, если форма «nikolskii» будет переописана с новой территории, то должны измениться ее название и диагностические признаки. Токсичность яда для мышей может быть использована в качестве нового диагностического признака.
Заключение
- Расчет точечных показателей токсикометрии (LD50 яда обыкновенных гадюк для мышей) по М. Л. Беленькому (1963) и с использованием современных версий пробит-регрессии не выявил существенных различий в конечных результатах.
- Использование обобщенных моделей, полученных методом максимального правдоподобия, позволяет изучать токсический процесс на всем диапазоне уровней воздействия, оценивать доверительные интервалы, анализировать различия в силе ядов на основе аппарата проверки статистических гипотез.
- Обобщенные линейные модели со смешанными эффектами (GLMEM) позволяют количественно оценить влияние комплекса сопутствующих случайных факторов на изучаемые процессы и выполнить прогноз изменения отклика в различных условиях воздействия.
- Токсичность яда может быть использована в качестве диагностического признака при переописании лесостепной формы обыкновенной гадюки.
Библиография
Бабич П. Н., Чубенко А. В., Лапач С. Н. Применение пробит-анализа в токсикологии и фармакологии с использованием программы Microsoft Excel для оценки фармакологической активности при альтернативной форме учета реакций // Современные проблемы токсикологии. 2003. № 4. С. 80–88.
Беленький М. Л. Элементы количественной оценки фармакологического эффекта . Л.: Медгиз, 1963. 152 с.
Мастицкий С. Э., Шитиков В. К. Статистический анализ и визуализация данных с помощью R . М.: ДМК Пресс, 2015. 496 с.
Международный кодекс зоологической номенклатуры. Издание четвертое . М.: Т-во науч. изданий КМК, 2004. 223 с.
Четвериков А. А. Линейные модели со смешанными эффектами в когнитивных исследованиях // Российский журнал когнитивной науки. 2015. Т. 2 (1). С. 41–51
Bakiev A. G., Böhme W., Joger U. Vipera (Pelias) berus nikolskii Vedmederya, Grubant und Rudaeva, 1986 – Waldsteppenotter // Handbuch der Reptilien und Amphibien Europas. Band 3/IIB: Schlangen (Serpentes) III. Viperidae. Wiebelsheim: AULA-Verlag, 2005. S. 293–309.
Bliss C. I. The calculation of the dosage-mortality curve // Annals of Applied Biology. 1935. Vol. 22. Issue 1. P. 134–167.
Burnham K. P., Anderson D. R. Model selection and multimodel inference: a practical information-theoretic approach. New York: Springer, 2002. 448 р.
Demidenko E. Mixed Models: Theory and Applications with R. 2nd ed. Wiley, 2013. 754 p.
Finney D. J. Probit Analysis. A statistical treatment of the sigmoid response curve. Cambridge: Cambridge University Press, 1947. 256 р.
Finney D. J. Probit Analysis. Cambridge: Cambridge University Press, 1971. 333 р.
Gelman A. Analysis of variance: Why it is more important than ever // The Annals of Statistics. 2005. Vol. 33. No 1. P. 1–53.
Milto K. D., Zinenko O. I. Distribution and Morphological Variability of Vipera berus in Eastern Europe // Herpetologia Petropolitana: Proceedings of the 12th Ordinary General Meeting of the Societas Europaea Herpetologica. St. Petersburg, 2005. P. 64–73.
Pinheiro J. C., Bates D. M. Mixed-Effects Models in S and S-PLUS. Springer, 2000. 530 p.
Robertson J. L., Savin N. E., Preisler H. K., Russell R. M. Bioassays with arthropods. Boca Raton: CRC Press, 2007. 224 p.
Zinenko O., Turcanu V., Stugariu A. Distribution and morphological variation of Vipera berus nikolskii Vedmederja, Grubant et Rudaeva, 1986 in Western Ukraine, The Republic of Moldova and Romania // Amphibia-Reptilia. 2010. Vol. 31. P. 51–67.
Благодарности
Авторы выражают благодарность Г. В. Еплановой, А. В. Павлову, М. В. Пестову, В. Г. Старкову и М. В. Ушакову за помощь в сборе образцов яда гадюк.
© 2011 - 2024