ОЦЕНКА ПАРАМЕТРОВ ВЫЖИВАЕМОСТИ ДЛЯ СОСТАВЛЯЮЩИХ НЕОДНОРОДНОЙ ПОПУЛЯЦИИ ПРОДУКТИВНЫХ ЖИВОТНЫХ: АНАЛИЗ ПРОБЛЕМЫ, ВАРИАНТЫ ПРИБЛИЖЁННОГО РЕШЕНИЯ Текст научной статьи по специальности «Математика»

МЕТОДЫ

УДК 636.2.034:591.139:519.654:51-33

ОЦЕНКА ПАРАМЕТРОВ ВЫЖИВАЕМОСТИ ДЛЯ СОСТАВЛЯЮЩИХ НЕОДНОРОДНОЙ ПОПУЛЯЦИИ ПРОДУКТИВНЫХ ЖИВОТНЫХ: АНАЛИЗ ПРОБЛЕМЫ, ВАРИАНТЫ ПРИБЛИЖЁННОГО РЕШЕНИЯ

1Черепанов Г.Г., 2Михальский А.И., 2Новосельцева Ж.А.

1ВНИИ физиологии, биохимии и питания животных, Боровск Калужской обл.;

2Институт проблем управления РАН, Москва, Российская Федерация

Цель данной работы - обоснование постановки задачи и предварительная отработка методик для оценки параметров выживаемости в популяции продуктивных животных, неоднородной по значениям средней и максимальной продолжительности хозяйственного использования в составляющих субпопуляциях или группах. Рассмотрены особенности двух методических подходов для такой оценки: 1) итерационные процедуры последовательного аналитического приближения по методу наименьших квадратов (МНК), 2) методы вычислительного эксперимента с использованием дискретной шкалы возраста. При оценке параметров выживаемости по данным возрастной структуры стада коров с использованием МНК необходимо учитывать влияние ряда «мешающих» факторов: 1) небольшое количество «точек» (пар значений интенсивности выбытия и номера лактации), соизмеримое с количеством определяемых параметров; увеличить количество точек можно, но для этого необходим переход к другой, более информативной системе регистрации и электронной архивации производственных данных; 2) статистическая надёжность оценки интенсивности выбытия выше на первых лактаци-ях, а на последующих она снижается, увеличивается разброс точек, нарушается требование равенства дисперсий отклонений, увеличивается риск получения «выскакивающих значений», вызывающих смещение оценок; 3) в период после исчерпания «короткоживущей» группы (когда интенсивность выбытия достигает единичного значения), возможно появление артефакт-ных значений её численности. При исследовании диаграмм, полученных с использованием «усечённого» распределения численности когорты по продолжительности продуктивной жизни (ППЖ) в режиме вычислительного эксперимента с дискретной шкалой возраста показано, что максимальная ППЖ евеличину средней ППЖПоказана возможность получения численных оценок ППЖ составляющих субпопуляций при условии, если значения с в показателе функции Гомпертца близки по величине, при этом разная длительность существования соответствующих когорт обусловлена, в основном, различиями в величине предэкспоненциально-го множителя. Согласно эмпирическим данным, эта величина коррелирует с уровнем жизнеспособности в ранней фазе репродуктивного периода. Использованные варианты решения поставленной задачи свидетельствуют о необходимости учитывать в исследованиях и в практической работе неоднородность популяций по параметрам выживаемости. Целесообразно также провести поиск дополнительной эмпирической базы, помимо динамики численности когорт, включая данные по частотам встречаемости генетических маркеров, ассоциированных с выживаемостью.

Ключевые слова: продолжительность жизни, неоднородные популяции, параметры выживаемости, методы исследования

Проблемы биологии продуктивных животных, 2017, 4: 81-95

Введение

оценки и ео направление ымтояние утраты жизнеспособности (Михальский и др., 1989; Михальский, 2002; Михальский, Цурко, 2016). «потенциалу долголетия», который для группы дочерей племенного быка можно оценить на первых лактациях (Черепанов и др., 2013; Черепанов, 2014). Оценка и учёт возможных эффектов гетерогенности стада по среднегруппо-вой и максимальной длительности хозяйственного использования могут иметь практическое значение, в том числе для совершенствования технологий разведения и селекции продуктивных животных.

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

Анализ проблемы. В течение длительного времени в демографических исследованиях считалось, что интенсивность смертности (количество умерших за некоторый отрезок времени, делённое на величину этого отрезка и число особей, доживших до начала этого отрезка) у человека и большинства животных, оцениваемая в группе одновременно родившихся генетически схожих особей, пропорциональна экспоненциальной функции от времени (возраста когорты, для человека приблизительно в период от 20 до 80 лет) с постоянным коэффициентом в показателе экспоненты. Это чисто стохастический механизм, из которого следует эмпирическое наблюдение, известное как экспоненциальная зависимость интенсивности смертности от возраста (Бенджамин Гомпертц, 1825). Хотя длительность хозяйственного использования продуктивных животных (продолжительность продуктивной жизни, length of productive life, longevity, herdlife) определяется их вынужденным выбытием из стада (выбраковкой по сумме причин), общие закономерности, характеризующие признак выживаемости, весьма схожи и для такого стада, и для популяций диких животных и человека. У племенных коров молочных пород при ежегодно проводимой бонитировке регистрируется продуктивность за 305 дней лактации, при анализе выживаемости возраст обычно выражают числом лактаций; при этом во многих зарубежных работах для анализа динамики выбытия используют функцию Вейбула (Ducrocq, Casella, 1996; Roxstrom et al., 2003), которая во многих практических ситуациях может быть аппроксимирована функцией Гомпертца.

В соответствии с этими классическими представлениями, появление любой патологии, приводящей к выбраковке коров, носит вероятностный характер, наступление очередной лактации рассматривается как случайное событие, в общем стаде формируются случайные потоки животных, переходящих из одного состояния в другое (Михайленко, 2015). Альтернативный вариант интерпретации — стадо неоднородно по выживаемости, т.е. состоит из нескольких групп с сильно различающимися значениями средней продолжительности хозяйственного использования и максимального регистрируемого в этой группе возраста (номера лактации). Согласно этой концепции, такие группы в каждом стаде представляют собой остатки когорт (животных одного года первой лактации) с разным начальным уровнем жизнеспособности, при этом внутри групп «работает» обычный вероятностный механизм выбытия. Под жизнеспособностью продуктивных животных здесь понимается мера функциональной мощности систем, обеспечивающих поддержание определённого уровня продуктивности, резистентности и репродуктивной способности на разных этапах жизни. То, что эта функциональная способность у всех животных снижается с возрастом, известно, а концептуальная новизна на данном этапе состоит в том, что в относительно молодом возрасте она, по всей вероятности, может быть представлена некой величиной — «потенциалом долголетия», которая определяет (в вероятностном смысле) длительность существования данной когорты (Черепанов и др., 2017).

Это два концептуально разных подхода — первый сосредоточен на оценке племенной ценности всех животных по одним и тем же критериям (безотносительно их потенциала жизнеспособности) и на борьбе с болезнями; второй подход нацелен на обеспечение максимально

возможного уровня «первичного здоровья» (Odent, 1986) и выявление группы «потенциальных долгожительниц». Различие в этих двух подходах имеет фундаментальный характер, поскольку здоровье и болезни — это в значительной мере разные сущности, понятие «здоровье» значительно шире, необходимость врачевания болезней можно рассматривать как неумение поддерживать здоровье в границах его возрастной нормы. Ещё И.П. Павлов утверждал, что медицина будущего — это профилактическая медицина, поэтому в XXI веке и ветврач, и селекционер должны уметь оценивать потенциал здоровья, повышать биологические возможности и степень надежности функциональных систем организма (особенно у потомства племенных малоплодных животных).

Этот аспект биомедицины в современном животноводстве становится всё более актуальным в связи с тем, что организм высокопродуктивных животных в критические периоды нередко (например, коровы с продуктивностью 12-14 тыс. кг в период раздоя) испытывает такие же перегрузки, которые у людей практикуются лишь в спорте высших достижений, но в этой области для формирования контингента рекордсменов проводятся длительные анализы и испытания и по их результатам «отсеиваются» сотни и тысячи претендентов, не подходящих по функциональным и конституциональным параметрам. Такое отсеивание не применяется для продуктивных животных, селекция молочного скота пока ведётся, в основном, по удою, статям и формам вымени, весь молодняк идёт «потоком» на воспроизводство стада, а по мере повышения его продуктивности (в процессе селекции и за счёт интенсивного кормления, применения биологически активных добавок и стимуляторов) по экспоненте снижается жизнеспособность, увеличивается выбраковка, страдает качество и безопасность получаемой продукции.

Рассмотрение в таком ракурсе проблем биологии продолжительности жизни в области животноводства означает перенос фокуса внимания с систем гомеостаза на системы гомеоре-зиса, проявляющегося в трендах изменений в организме на протяжении длительного времени. Трудность здесь в чрезвычайно высокой сложности и вариабельности объекта изучения; с этой проблемой сталкиваются все исследователи, работающие в области биологии продолжительности жизни, поскольку им по необходимости приходится иметь дело не с конкретными животными, а с популяциями и с длительными периодами времени «наблюдения», поэтому возникает необходимость прибегать к конструированию гипотез и концепций, для подтверждения (или опровержения) которых требуется организация длительных междисциплинарных исследований (Yashin, Iachine, 1997; Новосельцев, Новосельцева, 2007). В этой связи стоит указать, что аналогичный путь прошла генетика — 100 лет отделяет первую гипотезу о дискретной природе наследственности до расшифровки генетического кода и первых «работ руками» с генными конструкциями; от пионерской работы Гомпертца нас отделяет уже почти 200 лет, но успехи во многих направлениях биологии продолжительности жизни (в том числе в изучении популяционных структур выживаемости) весьма скромные, и это закономерно, поскольку объекты исследования в этой области имеют несравнимо большую степень системной сложности, в сравнении с объектами молекулярной генетики. Это отставание в большей степени сказывается на положении дел в области широкомасштабной селекции продуктивных животных, поскольку для человека проблема селекции не существует, а селекция лабораторных животных ведётся на немногочисленных популяциях чистых линий.

В настоящее время, при наличии общих сведений о средней по популяции продуктивности животных, отсутствуют принципы выявления генофондного стандарта породы, нет чёткого представления о таких понятиях, как потенциал продуктивности и резерв жизнеспособности, генетическая уникальность, недостаточно данных по оценке биоразнообразия для определения принадлежности особи к различным популяционным группам и т.д. Отсутствие фундаментальных знаний, стратегий и тактик сохранения генетических ресурсов животных и их рационального использования препятствует согласованным действиям на региональном и федеральном уровнях по управлению породообразовательными и продукционными процессами (Марзанов, 2007; Кузнецов, 2012, 2013; Черепанов, Марзанов, 2016).

За последние 25 лет численность поголовья дойных коров в целом в Российской Федерации сократилась в три раза (с 20 до 7 млн. голов), а средняя молочная продуктивность увеличилась приблизительно с 5 до 6 тыс. кг за лактацию. Серьёзным барьером на пути преодоления этого перекоса в современном промышленном животноводстве является антагонизм между уровнем продуктивности и жизнеспособностью маточного поголовья, выявляемый при анализе больших массивов эмпирических данных (Прошина, Лоскутов, 2011; Черепанов и др., 2017). Для преодоления этого антагонизма необходимо найти критерии для отбора родителей по показателям жизнеспособности и продуктивности потомства на возможно более ранней стадии, например, коров — по первым лактациям, так как успешность селекции, в особенности для малоплодных животных, сильно зависит от величины межгенерационного интервала.

Постановка задачи. Вопрос, возникающий в ряде практически важных приложений анализа продолжительности жизни, заключается в следующем. Если изучаемая популяция, в которой возрастная динамика интенсивности выбытия из популяции описывается функцией Гомпертца у(0 = В ехр(с^, неоднородна по выживаемости, то можно ли оценить численные значения параметров В и с для составляющих субпопуляций, а также их долю в общей популяции на основе общепопуляционных данных и, возможно, с привлечением дополнительной информации по другим источникам? Хотя аналитического решения этой задачи в такой общей постановке, по-видимому, не существует, но для получения оценок с приемлемой для практики точностью можно использовать методики расчётов, использующие итерационные процедуры последовательного приближения и оптимизации, а также методы вычислительного эксперимента.

Варианты приближённого решения задачи

Итерационная процедура получения оптимальных оценок

Исходные данные. В расчётах использовались эмпирические данные выбытия из стада коров холмогорской породы (Cherepanov et я1., 2013). В качестве исходной гипотезы предполагалось, что стадо состоит из пяти составляющих — остатков когорт с различными значениями параметра В при равенстве значений параметра с (обоснование этого условия приведено ниже). Задача — оценка начальных численностей групп и значений параметров В для составляющих стада, различающихся по средней продолжительности продуктивной жизни (ППЖ). Использовался нелинейный метод наименьших квадратов (МНК), алгоритм оптимизации Ньютона-Рафсона.

Методика анализа. Обозначим пропорции, с которыми пять субпопуляций представлены в популяции в начальный момент, и кумулятивный риск выбытия в каждой из субпопуляций через р,-(0) и Hj(t)=Bj/c*(exp(ct)-1) соответственно, где t — время, выраженное числом лакта-ций, ■ — последовательные номера лактаций. Суммарная функция дожития (текущая численность когорты) запишется в виде:

]=1

а интенсивность выбытия из общей популяции описывается выражением (2)

у(Т) = Е Р ('В ехр (ст) , (2),

]=1

где р>/(0 обозначает долю в общей популяции субпопуляции /; дожившей до момента времени I. Эта доля вычисляется по формуле

Р((,)= Р>(0)еХР(-И(<)) . (3)

£ р: (о)ехр (- И, (т))

,=1

Задача состоит в том, чтобы по эмпирическим данным и заданному значению параметра с найти оценки пропорций, с которыми пять субпопуляций представлены в обшей популяции в начальный момент. Для этого сначала необходимо оценить значения параметров В/, /=1,...5, используя соотношение (2) при Т, принимающем значения/ от 1 до 5 лактаций.

В качестве основы для такой оценки принимается теоретическое выражение для средней продолжительности жизни Т при заданных параметрах с и В в модели дожития Гомпертца (Кременцова, Горбунова, 2010):

1 Г в \г вЛ

Т = — ехр

с

* Е,

V с У

(4)

V с У

рда е у

где Е,(X) = I -йу - специальная функция (интегральная экспонента). Задавшись значе-

Л у

нием параметра с, для каждого Т/ определяется соответствующее значение параметра В/.

При найденных значениях В/ пропорции ^(0) находятся по эмпирическим данным нелинейным методом наименьших квадратов; при этом минимизируется сумма квадратов разностей между эмпирическими значениями интенсивности выбытия для последовательных лак-таций и значениями интенсивности выбытия, вычисленными по формуле (2) в те же моменты времени. Далее, при найденных значениях р/(0) вычисляется коэффициент детерминации Я2, который остаётся только функцией параметра с.

Максимизация коэффициента детерминации по с (с необходимым повторением описанной выше процедуры) даёт возможность получить соответствующее оптимальное значение параметров р/(0).

Заключение, которое можно сделать по результатам расчётов (с учётом сделанных предположений и в пределах разрешающей способности метода для данного набора данных), состоит в том, что в данном стаде отсутствуют представители короткоживущих групп со средней ППЖ 1 и 2 лактации, а также группа со средней ППЖ 4 лактации (рис. 1). Представители когорт со средней ППЖ 3 и 5 лактаций представлены практически в равных пропорциях; р(0) =51 и 49% соответственно. Коэффициент детерминации Я2 полученной зависимости для всего стада равен 0.835. Величина R2 показывает, какая часть (доля) дисперсии зависимой переменной объясняется рассматриваемой моделью. В случае линейной зависимости достаточно качественной считается модель с коэффициентом детерминации выше 0,8. Хотя для нелинейных моделей Я2 не является вполне адекватной характеристикой качества приближения рассматриваемой функции, в данном случае отклонение от линейной зависимости небольшое, поскольку эти отклонения проявляются более чётко по мере приближения к максимальной видовой продолжительности жизни, которая для коров в естественных условиях оценивается не менее чем 20 лактациями.

Рис. 1. Расчётные кривые возрастной зависимости интенсивности выбытия у(^) в неоднородной популяции коров холмогорской породы. Верхняя и нижняя линии - прогнозируемая методом наименьших квадратов интенсивность выбытия в двух субпопуляциях с указанными значениями Т. Чёрные точки - эмпирические значения интенсивности выбытия для стада численностью 1500 коров; показана аппроксимирующая их расчётная кривая для обшего стада, Я2= 0,835.

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

1) Небольшое количество «точек» (пар значений интенсивности выбытия и номера лактации) — как правило, не более 10; это накладывает определённые ограничения на количество оцениваемых параметров, так как при его увеличении свыше определённого уровня, специфического для каждого набора эмпирических данных, надёжность оценок снижается. В данном случае значения параметров Bs вычисляются с использованием условия, что все группы имеют одинаковые значения с. По эмпирическим данным оценивается 5 параметров (начальные доли групп в популяции). Оптимальное значение параметра с определяется из условия достижения максимума величины R2. Ввиду небольшого количества «точек», результаты проведенного расчёта следует рассматривать в качестве примера «работающего метода», а в дальнейшем необходимо проводить детальную его апробацию и исследовать возможности других методов. При этом будут иметь преимущество те методы, которые дополнительно дают значения стандартной ошибки оцененных параметров. Увеличить количество точек можно, но для этого необходим переход к более информативной системе регистрации и электронной архивации производственных данных (Черепанов и др., 2017).

2) В описанной процедуре минимизируется сумма квадратов разностей между прогнозируемыми и фактическими величинами по всему интервалу значений аргумента, при этом суммирование отклонений производится с равными весами для всех точек. Но в данном случае, ввиду того, что численность когорт с возрастом быстро снижается, а интенсивность выбытия вычисляется как относительная разность двух смежных значений, статистическая надёжность этой величины выше всего на первых лактациях, а на последующих она снижается, увеличивается разброс точек, нарушается требование равенства дисперсий отклонений. В связи с этим увеличивается риск получения «выскакивающих значений». Так, на рис. 1 критическое значение имеет значение ординаты для 5-й лактации, т.е. для остатка когорты коров, впервые отелившихся 5 лет назад. Если этот «выброс» связан с какими-то редкими форс-мажорными обстоятельствами, то такую точку обычно рекомендуют исключить и анализ проводить повторно.

3) Определённая проблема при использовании метода МНК для оценки параметров выживаемости коров возникает в случае существования «короткоживущих» субпопуляций, численность которых с возрастом быстро снижается, и если расчётная (при заданных значений В и с) интенсивность выбытия (-AS/S) по абсолютной величине становится больше 1, то прогнозируемая численность группы «короткожителей», т.е. функция дожития S(t), будет сначала небольшой отрицательной величиной (т.е. артефактной), с увеличением на последующих шагах и попеременным изменением знака. Следовательно, в применяемых алгоритмах необходимо предусматривать «обрыв» вычислений для короткоживущих когорт во избежание смещений в оценках параметров для общей популяции. Такие ситуации в демографических

работах не встречаются, так как в них рассматриваются возрастные интервалы в несколько десятков лет, но они возможны в популяциях продуктивных животных, если проводится «поперечный» анализ по данным возрастной структуры стада.

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

2) Вычислительный эксперимент (аппроксимация с использованием дискретных шкал и численных расчётов)

Предварительный анализ. При теоретическом анализе кривых выживаемости средняя продолжительность жизни Т определяется аналитической формулой (4) В этой формуле интегральная экспонента нужна для получения оценок с учётом предела видовой продолжительности жизни, но для многих практически важных ситуаций в расчётах можно использовать «усечённое» распределение численности когорты по продолжительности «жизни» — без учёта небольшого количества особей в области «хвоста» распределения.

Для проверки этого предположения использовали две методики расчётов — приближённое решение (аппроксимация с использованием дискретной шкалы возраста и численных расчётов) и там, где это возможно — точное аналитическое (с применением аналитических формул). Для получения исходного массива данных для приближённого решения была проведена серия из 98 реализаций расчета в формате Мюго80Й Ехсе11.

Исходные данные — значения В и с в функции Гомпертца у(0 = В ехр(с0 и начальная численность когорты N1. Для каждого значения В в диапазоне 0,04-0,3 использовали 7 значений с в диапазоне 0,2-0,8. Исследованный диапазон значений В соответствует эмпирическим данным по интенсивности выбытия коров из стада в период первой лактации, верхняя область диапазона значений с выбрана «с запасом» в связи с отсутствием достаточно многочисленных фактических данных по короткоживущим субпопуляциям (обычно регистрируется только возрастная структура дойного стада).

Использовались дискретные значения номеров лактации 7 = 1, 2, 3... 7тах. с соответствующими значениями возраста ^ и интенсивности выбытия из стада у, (y1=(N1-N2)/N1; у2=^2-N3)/N2 и т.д.) (табл. 1); — номер лактации; N1 — численность коров с номером лактации I;. п — количество выбывших на 7-ой лактации; п — сумма значений произведения п£ рав-

на средней продолжительности продуктивной жизни Т; ^ тах— последний номер лактации; в качестве аппроксимации используется «усечённое» распределение по продолжительности жизни в когорте без учёта количества особей в области «хвоста» распределения, не превышающего 1-5% от начальной численности.

Расчёт останавливали при достижении предельного значения ^тах, т.е. при появлении первого отрицательного значения величины (1 -у7) в случае короткоживущих субпопуляций или при достижении практически незначимой численности «долгожительниц». Значение средней продолжительности продуктивной жизни Т определяется суммированием п7%. по всем номерам лактации; под выбытием из стада понимается выбраковка «по сумме причин» (по причинам болезней, спада удоев молока или репродуктивной способности).

Выходные данные — значения 1тах , Т, у1, 1/у1 для данного сочетания параметров В и с — использовали для выявления взаимосвязей между параметрами выживаемости (табл. 2). Величину 1/у1 можно интерпретировать как среднегрупповой показатель жизнеспособности в ранней фазе репродуктивного периода (на первой лактации), т.е. как потенциал долголетия данной группы коров (Черепанов, 2014; Черепанов и др., 2016).

Таблица 1. Расчёт средней и максимальной продолжительности продуктивной жизни в когорте с параметрами В=0,14, с=0,3

г у 1-у N N041 п г*п

1 0,188 0,811 100 18,89 0,1889 0,1889

2 0,255 0,744 81,10 20,68 0,2068 0,4137

3 0,344 0,655 60,41 20,80 0,2080 0,6240

4 0,464 0,535 39,61 18,41 0,1841 0,7364

5 0,627 0,372 21,19 13,30 0,1330 0,6650

6 0,846 0,153 7,897 6,689 0,0668 0,4013

7 1,143 -0,143 1,209 1,381 0,0138 0,0967

8 1,543 -0,543 -0,173 -0,1731 -0,00276 -0,0213

9 2,083 -0,094 0,095 0,094 0,0009 3,1264

Т= Е г*п , = {3,126}

гтах 6 Т6= 3,11

Таблица 2. Фрагмент исходных данных для выявления взаимосвязей между параметрами выживаемости*

с В гтах Т * B/с У1 1/У1

0,2 0,2 7 2,93 0,586 1,000 0,24431 4,093

0,3 0,2 6 2,52 0,756 0,667 0,27002 3,703

0,4 0,2 5 2,22 0,888 0,500 0,29844 3,351

0,5 0,2 4 2,01 1,005 0,400 0,32985 3,032

0,6 0,2 4 1,85 1,11 0,333 0,36456 2,743

0,7 0,2 4 1,71 1,197 0,286 0,40293 2,482

0,8 0,2 3 1,56 1,248 0,250 0,44533 2,246

Примечание: * один из 14 циклов вычислений с варьированием параметров с и В.

Результаты и обсуждение. ри анализе ; величина е (рис. 2)

Рис. 2. Зависимость максимального возраста коров в когорте от средней продолжительности продуктивной жизни Т. По оси ординат - (номер последней лактации); по абсцисс - Т (среднее количество завершённых лактаций). Точки обозначают результаты 98 серий расчётов с варьированием численных значений параметров функции Гомпертца в интервале В = 0,04-0,3 и с = 0,2-0,8. Я2 — коэффициент достоверности аппроксимации (коэффициент детерминации).

Полученное соотношение можно использовать при поиске приближённого решения поставленной задачи и при выявлении взаимосвязей между параметрами выживаемости для составляющих субпопуляций с использованием пошаговых вычислений. Вычислительные эксперименты часто проводят при анализе сложных систем при неполных данных с целью получить ответ на вопросы такого типа: «Что будет, если...».

В данном случае использовали алгоритм последовательного поиска по следующей схеме: 1) задаётся несколько контрастных значений для когорт — составляющих субпопуляций; 2) по соотношению ~ 2определяются соответствующие величины Т; 3) путём скрининга всего массива данных, фрагмент которого приведен в табл. 2, по заданному значению Т определяются пары значений В и с (и соответствующие значения у1); 4) для каждой комбинации значений В и с после задания начальной численности когорт вычисляется «фракционная» функция Гомпертца (у) и последовательные величины численности от первой до последней лактации; 5) полученные значения численности когорт суммируются и по итоговым числен-ностям вычисляются значения функции Гомпертца для общей «модельной» неоднородной

популяции; 6) точечная диаграмма (зависимость интенсивности выбытия от возраста/номера лактации) итоговой «модельной» когорты сравнивается с аналогичной диаграммой для конкретного стада коров; 7) путём многократного «прогона» вычислений осуществляется поиск наиболее вероятной комбинации значений параметров.

Для примера 8-10 лактаций (групп) одна — «долгожительница», одна, относительно небольшая, быстро выбывает из стада (^тах ~ 2) , и две (^тах ~ 4-8) необходимо иметь для того, чтобы «срединных»и общей формы

Для оценки точности алгоритма приближённого расчёта были сопоставлены две функции, вычисленные по результатам численного эксперимента и по аналитическим формулам. Вычислялось произведение параметра с*Т в зависимости от отношения В/с , а также значение Т в зависимости от величины 1/у1. При расчёте Т учитывалось, что речь идёт о продолжительности продуктивной жизни, которая короче видовой продолжительности жизни и вычисляется по формуле (5)

Т = ^ехр | — !* с I с

Е' — ) Е ( — еХР ^

(5)

где ^тах - максимальная длительность продуктивной жизни, принятая для четырёх субпопуляций равной 2, 5, 8 и 10 лактаций.

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

с*Т =/(В/с)

Аппроксимация

В/с

3

2.5 2

1.5 1

0.5

0

с*Т

Аналитическое решение с *Т=ехр(В/с)*Е(В/с)

теоретическая зависимость аппроксимация

0 0.2 0.4 0.6 0.8 1 1 2 1.4 1.6 В/С

Рис. 3а. Зависимость произведения с*Т от величины отношения В/с. Слева - результаты численного эксперимента (аппроксимация), справа - прогноз по аналитической формуле (4). В численном эксперименте использована дискретная шкала последовательных номеров лактаций и «усечённое» распределение численности когорты по продолжительности жизни (без учёта небольшого количества особей в области «хвоста» распределения). Точки -98 пар значений В и с.

Аппроксимация

T=f(1/yi)

1/yi

Аналитическое решение

T= с1 exp(e'c/(c/y1))*[Ei(e'c/(c/y1))-Ei(e("n'^'1c/(c/y1)] T

8 7 6 5 4 3 2 1

0 2 4 6 8 10 12 14 16 18 20 22

1/y i

♦ c= 0.2 O c= 0.4

• c= 0.6 ¿ c = 0.8

Рис. 3б. Зависимость средней продолжительности продуктивной жизни Т от величины 1/У\ (у1 - интенсивность выбытия на первой лактации). Слева - результаты численного эксперимента (аппроксимация), справа - прогноз по аналитической формуле (линии) в сопоставлении с результатами численного эксперимента (значки).

Численный пример анализа популяции, неоднородной по параметрам выживаемости.

Апробация методики проведена на том же массиве эмпирических данных, как и в разделе (1). Задачей было оценить значения параметров В, с и начальной относительной численности четырех субпопуляций (групп), с тем, чтобы прогнозируемая функция у = ДО для суммарной модельной популяции была близка к эмпирическим данным. Исходно задавались значения Т = 1, 2, 3, 4 лактации для четырёх групп; процедуры оценивания параметров этих групп проводились последовательно по убывающим значениям ^тах = 2Т в направлении с «конца к началу», поскольку при этом предыдущие когорты, по приведенным выше оценкам, исчерпывались до достижения данного значения Т. По заданным значениям Т путём скрининга всего массива данных, фрагмент которого приведен в табл. 2, определялись пары значений В и с и «методом проб и ошибок» делались попытки восстановить эмпирическую кривую интенсивности выбытия в общей популяции.

На этом этапе выявилась закономерность, имеющая существенное значение для всего последующего анализа. Оказалось, что такое приближение возможно при условии, если значения с для субпопуляций, имеющих значения Т = 2, 3 и 4, близки по величине к одному значению — 0,2 (как и для общей популяции), так что различная длительность существования соответствующих когорт обусловлена, в основном, различиями в величине В. Ранее на большом массиве эмпирических данных было показано, что различия коров по величине ППЖ в значительной степени обусловлены разным уровнем их конститутивной резистентности в ранней фазе репродуктивного периода (Черепанов, 2014). С другой стороны, схожие величины параметра с для разных групп, по-видимому, могут быть отражением общих породных особенностей, географического фона и условий эксплуатации молочного стада.

В проведенном численном эксперименте противоположно направленные вариации в значениях с для короткоживущей группы не оказывали большого влияния на степень приближения результирующей кривой к фактическим данным, если начальная численность этой когорты составляла небольшую долю от соответствующей величины в общей когорте (в данном случае — 5-6%) (рис. 4). Поэтому можно предположить, что нулевое значение начальной доли «короткоживущей» составляющей в первом варианте расчётов (рис. 1) получилось в силу того, что небольшая величина этого параметра оказалась вне предела «разрешающей способности» метода для данного набора эмпирических данных. В проведенном эмпирическом подборе кривых для «субпопуляций» учитывалась выше отмеченная особенность — более высокая статистическая надёжность оценок на первых лактациях, поэтому результирующая кривая проходит более близко к эмпирическим значениям, чем на рис. 1 .

Рис. 4. Два варианта расчётных кривых возрастной зависимости интенсивности выбытия уф в неоднородной популяции коров холмогорской породы (по данным вычислительного эксперимента). По оси ординат - интенсивность выбытия уф; по оси абсцисс - номера лактаций. Тонкие сплошные линии - прогнозируемые функции у для четырёх субпопуляций; о - эмпирические данные; х - расчётные значения для общей популяции; жирная линия -прогнозируемая функция для общей популяции, упрогноз.

На данном этапе не ставилось целью получение точных значений параметров с оценкой статистической погрешности; немаловажно здесь то, что оба использованных варианта решения предсказывают наличие в стаде нескольких достаточно крупных групп (две в первом варианте и 3 во втором, с суммарной долей в стаде не менее 90%), различающихся по величине средней и максимальной продолжительности продуктивной жизни. Что собой представляют такие группы, и какие факторы служат причиной появления неоднородности популяций по параметрам выживаемости — это вопросы для последующих детальных исследований.

Общее обсуждение и выводы

Два методических приёма, рассмотренные в разделах 1) и 2), как было отмечено выше, имеют свои достоинства и недостатки. Во второй методике значения показателя В и процентной доли для составляющих субпопуляций имеют ориентировочный характер, поскольку не проведена процедура оптимизации численных оценок, а в первом разделе требуют доработки вопросы, связанные с неоднородностью дисперсий отклонений, влиянием «выскакивающих» значений и возможным появлением артефактных величин после исчезновения короткоживу-щих групп. Тем не менее, оба варианта решения дали дополнительную аргументацию о возможности и необходимости учитывать в научных исследованиях и в практической работе неоднородность популяций по выживаемости, т.е. существование групп, различающихся по величине средней и максимальной продолжительности хозяйственного использования. Поэтому в дальнейшей работе в этом направлении следует развивать комплексный подход с использованием различных методических приёмов.

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

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

В качестве численного примера, рассмотрим ситуацию с двумя аллелями — А и В (необязательно ассоциированными с выживаемостью), которые выявляются в начальном возрасте tj с частотами 0,7 и 0,3 соответственно в выборке объёмом 1000 (т.е. число носителей А и В — 700 и 300 соответственно). Если составляющие выборку когорты однородны по выживаемости (т.е. имеют одинаковые значения параметров функции, описывающей динамику интенсивности выбытия), количество оставшихся носителей А и В в возрасте t2 в среднем снизится в одинаковой пропорции, и частоты встречаемости не изменятся. Если же когорта неоднородна, например, носители А — долгоживущие, а носители В — короткоживущие, количество оставшихся будет разным, например, если количество носителей А снизится в 2 раза, а В — в 4 раза, то останется 350 и 75 особей соответственно. В случайной выборке из этой оставшейся смешанной группы соответствующие частоты будут равны 350/425=0,82 и 75/425=0,18, т.е частота А увеличится, а частота В снизится, в сравнении с начальными значениями в момент tj. В данной ситуации это связано не с возрастным полиморфизмом, как может показаться на первый взгляд, а с разной скоростью выбытия носителей аллелей А и В в двух когортах — составляющих общей неоднородной популяции.

Таким образом, если в случайно взятой выборке коров данного возраста обнаруживается сдвиг частот встречаемости генетического маркера, в сравнении с группой первотёлок, то это можно считать вероятным эмпирическим свидетельством неоднородности стада по параметрам выживаемости (средней и максимальной III1Ж), — например, одна когорта исчезает через 2 года, а вторая доживает до 6-8. Эмпирические данные по возрастным сдвигам частот встречаемости аллелей и генотипов в популяциях коров имеются, хотя такие работы немногочисленны (Иванов и др., 2017).

В целом, в данной работе проанализированы возможные варианты приближённой оценки параметров выживаемости в неоднородной популяции продуктивных животных, т.е. состоящей из нескольких субпопуляций или групп, различающихся по значениям средней и максимальной продолжительности хозяйственного использования. Показано, что на первом этапе работы в этом направлении можно получить приближённые оценки параметров В и с в функции Гомпертца y(t)=B exp(ct) по неполным данным — по эмпирической функции интенсивности выбытии для общей популяции с использованием алгоритмов, минимизирующих расхождения между значениями прогнозируемой функции для суммарной модельной популяции и эмпирическими данными.

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

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

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

ЛИТЕРАТУРА

1. Иванов В.А., Марзанов Н.С., Елисеева Л.И., Таджиев К.П., Марзанова С.Н. Генотипы пород крупного рогатого скота и качество молока // Проблемы биологии продуктивных животных. - 2017. - № 3. - С. 48-65.

2. Кременцова А.В., Горбунова Н.В. Роль окружающей среды в динамике распределения продолжительности жизни // Автоматика и телемеханика. - 2010. - № 8. - С. 121-133.

3. Кузнецов В.М. Племенная оценка животных: прошлое, настоящее, будущее (обзор) // Проблемы биологии продуктивных животных. - 2012. - №4. - С. 18-57.

4. Кузнецов В.М. Разведение по линиям и голштинизация: методы оценки, состояние и перспективы // Проблемы биологии продуктивных животных. - 2013. - № 3. - С. 25-79.

5. Марзанов Н.С. Генетические ресурсы животных - это национальная ответственность // Известия ТСХА. - 2007. - Вып. 5. - С. 135-141.

6. Михальский А.И., Петровский А.М., Яшин А.И. Теория оценивания неоднородных популяций. - М.: Наука, 1989.

7. Михальский А.И. Методы анализа гетерогенных структур в популяции. - М. : ИПУ РАН институт проблем управления, 2002. - 62 с.

8. Михальский А.И., Цурко В.В. Анализ рисков с учётом гетерогенности и сочетанной заболеваемости // Мат. девятой межд. конф. «Управление развитием крупномасштабных систем MLSD'2016». - М.: ИПУ РАН, 2016. -Том 2. - C. 388-390.

9. Новосельцев В.Н., Новосельцева Ж.А. Моделирование истории жизни и баланса ресурсов // В кн.: Геронтология in silico: становление новой дисциплины (Ред. Г.И. Марчук, В.Н. Анисимов, А.А. Романюха, А.И. Яшин). - М.: БИНОМ. Лаборатория знаний, 2007. - С. 148-174.

10. Прошина О., Лоскутов Н. Воспроизводство стада: потерянная страница // Животноводство России. -2011. - № 9. - 40-41.

11. Черепанов Г.Г. Обоснование концепции о ключевой роли конститутивной резистентности для жизнеспособности и длительности использования высокопродуктивных животных // Проблемы биологии продуктивных животных. - 2014. - № 4. - С. 5-34.

12. Черепанов Г.Г., Михальский А.И. Анализ возможных подходов для оценки жизнеспособности и прогнозирования долголетия у высокопродуктивных животных // Проблемы биологии продуктивных животных. - 2016. - № 1. - С. 5-25.

13. Черепанов Г.Г., Харитонов Е.Л., Макар З.Н., Михальский А.И., Новосельцева Ж.А. Анализ возможных подходов для преодоления антагонизма между уровнем продуктивности и жизнеспособностью маточного поголовья при использовании интенсивных технологий // Проблемы биологии продуктивных животных. - 2017. - № 1. - С. 5-27.

14. Черепанов Г.Г. Марзанов Н.С. Пилотные модели продукционных процессов для использования в крупномасштабных системах управления в животноводстве // В сб.: Материалы девятой международной конференции: «Управление развитием крупномасштабных систем (MLSD'2016). - М.: ИПУ РАН, 2016. - Т. 2. - С. 411-413.

15. Cherepanov G.G., Bogdanova N.A., Makar Z.N. Spécial traits of the age dynamics of milk production in dairy cows in relation to their viability // Russian Agricultural Sciences. - 2013. - Vol. 39. - No. 5-6. -Р. 442-445.

16. Ducrocq V., Casella G. A Bayesian analysis of mixed survival models // Genet. Sel. Evol. - 1996. - Vol. 28. -Р. 505-529.

17. Roxstrom A., Ducrocq V., Strandberg E. Survival analysis of longevity in dairy cattle on a lactation basis // Genet. Sel. Evol. - 2003. -Vol. 35. - P. 305-318.

REFERENCES

1. Cherepanov G.G., Bogdanova N.A., Makar Z.N. Special traits of the age dynamics of milk production in dairy cows in relation to their viability. Russian Agricultural Sciences, 2013, 39(5-6): 442-445.

2. Cherepanov G.G. [Substantiation of the concept of the key role of constitutive resistance for the viability and duration of use of highly productive animals]. Problemy biologiiproductivnykh zhivotnykh - Problems of Productive Animal Biology. 2014, 4: 5-34.

3. Cherepanov G.G., Kharitonov E.L., Makar Z.N., Mikhal'skii A.I., Novosel'tseva Zh.A. [An analysis of possible approaches to overcome the antagonism between the level of productivity and the viability of the breeding stock under using intensive technologies]. Problemy biologii productivnykh zhivotnykh - Problems of Productive Animal Biology. 2017, 1: 5-27.

4. Cherepanov G.G., Mikhal'skii A.I. [Analysis of possible approaches for assessing the viability and prediction of longevity in highly productive animals]. Problemy biologii productivnykh zhivotnykh - Problems of Productive Animal Biology. 2016, 1: 5-25.

5. Cherepanov G.G., Marzanov N.S. [Pilot models of production processes for use in large-scale management systems in animal husbandry]. In: Mat. devyatoi mezhd. konf. « Upravlenie razvitiem krupnomasshtabnykh sistemMLSD'2016» (Mat. 9th Int. Conf.: Management of large-scale systems development - MLSD'2016). Moscow: IPU RAS Publ., 2016, 2, P. 411-413.

6. Ducrocq V., Casella G. A Bayesian analysis of mixed survival models. Genet. Sel. Evol. 1996, 28: 505-529.

7. Ivanov V.A., Marzanov N.S., Eliseeva L.I., Tadzhiev K.P., Marzanova S.N. [Genotypes of breeds of cattle and quality of milk]. Problemy biologii productivnykh zhivotnykh - Problems of Productive Animal Biology. 2017, 3: 48-65.

8. Krementsova A.V., Gorbunova N.V. [The role of the environment in the dynamics of the distribution of life expectancy]. Avtomatika i telemekhanika - Automation and telemechanics. 2010, 8: 121-133.

9. Kuznetsov V.M. [Breeding animals rating: past, present, future: a review]. Problemy biologii productivnykh zhivotnykh - Problems of Productive Animal Biology. 2012, 4: 18-57.

10. Kuznetsov V.M. [Line breeding and golshtinization: assessment methods, current state and prospects]. Problemy biologii productivnykh zhivotnykh - Problems of Productive Animal Biology. 2013, 3: 25-79.

11. Marzanov N.S. [Animal genetic resources are national responsibility]. Izvestiya TSKhA - Proc. TSKhA. 2007, 5: 135-141.

12. .Mikhal'skii A.I., Petrovskii A.I., Yashin A.I. Teoriya otsenivaniya neodnorodnykh populyatsii (Theory of estimating non-homogeneous populations). Moscow: Nauka Publ., 1989, 128 p.

13. Mikhal'skii A.I. Metody analiza geterogennykh struktur vpopulyatsii (Methods for analyzing heterogeneous structures in a population). Moscow: IPU RAS Publ., 2002, 62 p.

14. Mikhal'skii A.I., Tsurko V.V. [Risk analysis taking into account heterogeneity and combined morbidity]. In: Mat. devyatoi mezhd. konf. «Upravlenie razvitiem krupnomasshtabnykh sistem MLSD'2016» (Mat. 9th Int. Conf.: Management of large-scale systems development MLSD'2016). Moscow: IPU RAS Publ., 2016, 2: 388-390.

15. Novosel'tsev V.N., Novosel'tseva Zh.A. [Modeling life history and balance of resources]. In: Gerontologiya in silico: stanovlenie novoi distsipliny (Gerontology in silico: the emergence of a new discipline. Eds G.I. Marchuk, V.N. Anisimov, A.A. Romanyukha, A.I. Yashin). Moscow: BINOM Laboratoriya znanii Publ., 2007, p. 148-174.

16. Proshina O., Loskutov N. [Reproduction of the herd: lost page]. Zhivotnovodstvo Rossii - Animal Husbandry in Russia. 2011, 9: 40-41.

17. Roxstrom A., Ducrocq V., Strandberg E. Survival analysis of longevity in dairy cattle on a lactation basis // Genet. Sel. Evol. - 2003. -Vol. 35. - P. 305-318.

Evaluating parameters of survival for the components of nonhomogeneous population of productive animals: analysis of problem, development of approximate solution

:Cherepanov G.G., 2Mikhalskii A.I., 2Novoseltseva Zh.A.

'Institute of Animal Physiology, Biochemistry and Nutrition, Borovsk, Russian Federation; 2Institute of Control Sciences RAS, Moscow, Russian Federation

ABSTRACT. In recent years, in the field of lifetime biology, the concept of the need to take into account the heterogeneity of populations in terms of survival parameters has been developed. The purpose of this work is to justify the formulation of the problem and preliminary development of techniques for the approximate evaluation of survival parameters in a population of productive animals that is inhomogeneous in terms of mean and maximum length of productive life in subpopulations or groups. Features of two methodological approaches for such estimation are considered: 1) iterative procedures of sequential analytical approximation by the least squares method (LSM), 2) methods of computational experiment using the discrete age scale. When assessing survival parameters from the age structure of a herd of cows using LSM, it is necessary to take into account the influence of a number of "disturbing" factors: 1) a small number of "points" (pairs of values of culling intensity and lactation number) commensurable with the number of parameters to be determined; the increase in the number of points is possible, but this requires a transition to another, more informative system of registration and electronic archiving of production data; 2) the statistical reliability of estimating the culling intensity is higher on the first lactations, and on subsequent ones it decreases, the scatter of points increases, the requirement of equality of deviation variances is violated, the risk of obtaining "pop-up values" causing increase in estimates bias; 3) in the period after the exhaustion of the "short-lived" group (when the predicted culling intensity reaches a single value), the appearance of artifact values of group size is possible. In the study of the diagrams obtained using the "truncated" distribution of the cohort population by the length of productive life (LPL) in the mode of computational experiment with a discrete age scale, it was shown that the maximum LPL in the model subpopulations is approximately 2 times larger than the mean LPL. The possibility of obtaining numerical estimates of LPL for subpopulations is shown under the condition that the values of c in the exponent of the Gompertz function are close in size, while the different lifetimes of the corresponding cohorts are mainly due to differences in the magnitude of the pre-exponential factor. According to empirical data, the value of this factor correlates with the level of constitutive resistance in the early phase of reproductive period. The variants used to solve this problem testify the need to take into account the inhomogeneity of populations in survival parameters in studies and in practical work. It is also advisable to search for ways of obtaining an additional empirical base, in addition to the dynamics of the cohort number, including data on the frequency of occurrence of genetic markers associated with survival.

Keywords: cows of dairy breeds, length ofproductive life, heterogeneous populations, survival parameters, research methods

Problemy biologii productivnykh zhivotnykh - Problems of Productive Animal Biology, 2017, 4: 81-95

Поступило в редакцию: 23.07.2017 Получено после доработки: 03.09.2017

Черепанов Геннадий Георгиевич, д.б.н., г.н.с., т. 8(905)642-03-99; [email protected];

Михальский Анатолий Иванович, д.б.н., к.т.н., т. (915)199-55-26; [email protected], [email protected];

Новосельцева Жанна Анатольевна, к.т.н., с.н.с., т. (495)334-88-91; [email protected]