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

Исследователи давно признали, что математическое моделирование и компьютерные технологии могут помочь клиницистам в принятии решений, помогая им осуществлять процесс «просеивания» возможных диагнозов и симптомов с целью выработки тактики лечения и/или прогнозирования исхода болезни [1]. Первые попытки анализа биомедицинской информации с использованием компьютерных технологий были предприняты в 60-х годах прошлого столетия. Однако сегодня в медицинской практике используются лишь немногие компьютерные прогностические модели, многие из которых неточны и обладают слабой клинической достоверностью, не обладают универсальностью, а их эффективность низка [2, 3].

Одной из немногих весьма популярных среди онкологов программ для таких прогнозов является бесплатный Интернет-ресурс Adjuvant! Online [https: //www.adjuvantonline.com/index.jsp]. В этой программе для прогнозирования рисков используется метод the Prognostic Factor Impact Calculator, основанный также на методологии байесовских сетей (БС) [4], но в описании методологии вычисления нет, поэтому проверить надежность работы the Prognostic Factor Impact Calculator, например узнать величину его AUC, к сожалению, невозможно. Также к числу недостатков этого ресурса следует отнести тот факт, что количество параметров пациентов, на основании которых делается прогноз, очень мало — всего 7 наименований, и ввести дополнительные параметры невозможно, так как программа не имеет соответствующей опции. Кроме того, последняя версия программы написана в 2005 г. и добавить новые данные о пациентах, например информацию об экспрессии рецептора HER-2/neu, в эту программу невозможно. Этот Интернет-ресурс работает в бо́льшей степени как «black-box», в результате чего довольно трудно понять как качество и надежность прогнозов, так и обоснованность выводов.

Тем не менее вероятностные модели на основе БС становятся все более популярными в медицине [5—9]. Аппарат Б.С. представляет собой удачную комбинацию теории вероятности и теории графов, позволяющую объединять различные клинические стратегии в лечении больных и различные типы таких данных, как клинико-патоморфологические и огромные базы генетического анализа опухолей [10].

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

Список параметров, применяемых в онкологии, постоянно пополняется, происходит поиск новых молекулярно-биологических прогностических и предиктивных факторов.

Так, X. Jiang и соавт. [2] описали метод EBMC—S, основанный на технологии БС, и его применение для прогнозирования выживаемости больных раком молочной железы (РМЖ). Авторы используют структурную оптимизацию байесовской сети с использованием некоторого приближенного критерия оптимизации — скора. Эффективность предложенного метода авторы показывают, используя базу данных больных РМЖ, и сравнивая свои прогнозы по выживаемости с результатами стандартного метода пропорциональных рисков Кокса [11] и метода RSF (the random survival forest method) [12]. Показано, что метод EBMC—S существенно лучше метода пропорциональных рисков Кокса и несколько лучше метода RSF. Тем не менее у этого метода выявили два недостатка. Во-первых, величина AUC (площадь под ROC-кривой) весьма невелика и составляет величины для прогнозирования на разных временных интервалах от 0,61 для одного года до 0,77 для 15 лет. Во-вторых, можно обнаружить результат (рис. 1), который противоречит здравому смыслу — с ростом времени качество прогнозирования выживаемости увеличивается. В этой связи необходима дальнейшая работа по улучшению качества прогнозирования рисков исхода для пациенток с РМЖ.

Сегодня стратегия и тактика лечения пациентов, прогноз для жизни больных раком осуществляются на основании так называемого «молекулярного портрета опухоли». В качестве нового прогностического маркера мы предложили и исследовали белок YB-1.

Белок YB-1 принадлежит к классу ДНК/РНК-связывающих белков с эволюционно-консервативным доменом холодового шока [13]. Белок YB-1 тесно связан с опухолевой прогрессией: он участвует в возникновении множественной лекарственной устойчивости [14], активирует пролиферацию и миграцию клеток [15, 16], связан с PI3K/Akt/mTOR-сигнальным путем [17, 18].

Исходя из роли белка YB-1 в клетке, исследователи предположили, что его повышенная экспрессия и/или ядерная локализация может быть показателем агрессивности течения различных онкологических заболеваний. И действительно, было показано, что ядерная локализация белка YB-1 коррелирует с худшим прогнозом больных немелкоклеточным раком легкого и раком яичников [19, 20]. Нами ранее показано, что повышенная экспрессия мРНК гена YB-1, исследуемая методом ПЦР в реальном времени и методом полуколичественного ОТ-ПЦР коррелирует с увеличением частоты метастазирования РМЖ и худшей безрецидивной выживаемостью больных [21]. Таким образом, можно заключить, что использование показателей уровня экспрессии гена YB-1 и его локализации могут существенно влиять на прогноз исходов для больных РМЖ и усилить прогностическую способность создаваемой экспертной системы персонифицированной медицины.

Цель исследования — применение технологии БС для прогнозирования риска прогрессирования, а также общей и безрецидивной выживаемости больных РМЖ, используя различные прогностические параметры, в том числе показатели экспрессии и локализации YB-1.

На основе построенной и оптимизированной БС в данной работе были определены критические параметры для предсказания риска прогрессирования и смерти. Экспрессия мРНК гена YB-1 в опухолевой ткани вошла в число критических параметров, определяющих риск наступления летального исхода в первые годы после проведения хирургического лечения. Построенные на основе обученных БС гистограммы риска, прогнозирующие индивидуальную вероятность исхода заболевания, могут помочь врачу-онкологу при принятии решения по тактике лечения больных РМЖ.

В ретроспективное исследование вошли пациентки с верифицированным диагнозом РМЖ I—IV стадии в возрасте от 30 до 84 лет. Больные проходили лечение в Онкологическом центре ЦКБ им. Н.А. Семашко ОАО «РЖД» в период с 1992 по 2008 г. База данных включала стандартные параметры (варианты VAR001-VAR0051) 323 больных, необходимые для принятия решения о противоопухолевом лечении РМЖ (табл. 1). Некоторые параметры, такие, например, как степень дифференцировки опухоли и др. (см. табл. 1), оказались частично недоступны и поэтому каждый столбец «Процент заполнения данных» в сумме может не составлять 100%. Последний мониторинг состояния пациенток мы осуществляли в 2013 г.

Экспрессию мРНК гена YB-1 (метод ОТ-ПЦР) и локализацию белка YB-1 (иммуногистохимический метод) определяли в опухолевом материале, полученном при мастэктомии по методике, описанной ранее [21].

Для описания параметров состояния пациенток были определены следующие значения: состояние без прогрессирования, прогрессирование заболевания, метастатическая болезнь (ТхNхМ1) и смерть. На основе известных состояний мы исследовали две конечные точки (КТ): КТ1-П — прогнозирование прогрессирования заболевания и КТ2-С — прогнозирование смерти пациентки. Для каждой КТ мы определили два исхода: КТ1-П принимает значение «0», если заболевание не прогрессирует, и значение «1», если заболевание прогрессирует, КТ2-С принимает значение «0», если пациентка жива, и значение «1», если пациентка умерла. Соответствие состояний пациенток и значений КТ представлены в табл. 2.

Нами было предложено прогнозировать КТ с момента проведения у пациентов операции. Так как время проведения операции было разным, а время мониторинга за их состоянием одинаковое, то необходимо было провести стратификацию пациенток по временным интервалам. Для этого в базу данных ввели новый параметр «status praesens» (SP), который принимал значение состояния пациентки на момент заданного временного срока. Значение этого параметра определяли с учетом информации о состоянии пациентки при ближайшем мониторинге по правилам, указанным в табл. 3.

Было исследовано соотношение значений КТ1-П и КТ2-С для временных сроков от 1 года до 6 лет. В табл. 4 представлено число пациенток по годам для соответствующих значений К.Т. Из этой таблицы видно, что общее число пациенток для выбранных КТ существенно падает для 5 и 6 лет по сравнению с интервалом 4 года, и поэтому для дальнейшего анализа мы выбрали временной интервал 4 года, несмотря на то что в онкологии принятым временным интервалом является 5-летний срок.

Для дальнейшей работы из базы данных было выбрано 32 параметра (см. табл. 1): основные клинические и молекулярно-биологические данные, информация о видах лечения. Непрерывные параметры были разбиты на интервалы, и проведена их дискретизация. Для этого области значений мы разбили на несколько интервалов, так что все переменные в базе данных принимали целочисленные значения. На основе данных параметров для каждого из двух КТ были построены БС с наивной топологией (см. ниже). Сети отличались только корневыми узлами, определяющими значения КТ.

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

Выбор топологии БС существенно влияет на значения условных вероятностей. В настоящей работе мы использовали наивную топологию сети (рис. 2), т. е. у всех узлов есть один общий родитель (корневой узел «A»). В случае персонифицированной медицины в качестве корневого узла «A» используется КТ.

Данная топология имеет ряд преимуществ. Во-первых, в данной модели получается простое выражение для полной вероятности. Во-вторых, достигается наименьшее количество таблиц условных вероятностей. В-третьих, отсутствует влияние переменных друг на друга, что уменьшает минимально необходимое количество данных для обучения. Это особенно важно при работе с малыми базами данных пациентов. Построение наивных БС осуществлялось разработанной нами программой NetGen (Network Generating).

Для сравнения прогнозируемой способности БС мы применяли метод ROC-кривых [22—23]. Для численного сравнения качества прогноза БС на исследуемых базах данных использовалось значение площади под ROC-кривой (AUC): чем ближе значение AUC к 1, тем выше качество прогноза обученной сети. Для построения ROC-кривой применяли метод «исключением по одному» (рис. 3). База данных была представлена в виде таблицы, в ячейках которой стоят целые числа, а каждая строка является данными для определенного пациента. При работе алгоритмов обучения и опроса строки с пропущенным значением КТ игнорировали. На каждом шаге построения ROC-кривой из базы данных вычеркивали одного пациента, на оставшихся пациентах происходило обучение БС, и далее вычеркнутого пациента опрашивали на обученной БС, т. е. по параметрам вычеркнутого пациента мы прогнозировали исход (КТ) и определяли соответствующую условную вероятность. По окончании каждого шага формировали ROC-таблицу, содержащую в каждой строке реальный исход и условную вероятность наступления этого исхода. Далее таблицу сортировали по условной вероятности, и на ее основе строили ROC-кривую и вычисляли значение AUC.

Поскольку мы использовали наивную топологию, алгоритмы обучения и опроса БС были проще, чем аналогичные алгоритмы для общего случая [24—26]. Алгоритм обучения представлял собой прямое использование определения условной вероятности. Данный алгоритм был точным и не имел случайных факторов, т. е. получался один и тот же результат для любого запуска на определенной сети и заданной базе данных. Все таблицы условных вероятностей были представлены в форме потенциалов. Использование наивной топологии заметно упрощает алгоритм опроса до перемножения всех потенциалов с последующей нормировкой на единицу.

Оценку качества прогноза методом вычисления AUC «исключением по одному» осуществляли с помощью разработанной нами программы ANN (AUC for Native Network).

Одной из проблем, возникающих при использовании БС, является выделение множества переменных, имеющих критическое значение для прогнозируемой способности сети, и устранения переменных, слабо или отрицательно влияющих на нее. Пусть в качестве численной оценки прогнозируемой способности заданной сети при обучении на конкретной базе выступает некоторый критерий f (X) (в данном случае X — множество переменных рассматриваемой сети). В нашем случае роль f (X) играет величина AUC. Тогда оптимизация сети заключается в поиске множества переменных, которые, будучи использованы для построения БС при работе с заданной базой данных, дадут максимально большое значение f (X). В нашем случае оптимальной будет та БС, которая даст максимальное значение AUC.

Предлагаемый алгоритм оптимизации сети основан на анализе возможных сетей, составленных из различных комбинаций переменных исходной сети. В каждую такую комбинацию обязательно входит корневая переменная исследуемой БС, соответствующая заданной К.Т. Чтобы исключить влияние топологии на величину f (X) и снизить время расчетов, рассматриваются лишь сети с наивной топологией. В таких сетях корневая переменная сети является родительской для всех остальных переменных, эти переменные не связаны между собой. Блок-схема алгоритма выявления искомого множества переменных представлена на рис. 4.

На рис. 4 видно, что алгоритм является итерационным, на каждой итерации вычисляется значение f (X0) исходной сети, включающей все переменные (узлы) этой сети, а также значения f (X) всех сетей, получающихся из исходной сети исключением одной переменной, кроме корневой. Из этих сетей выбирается та сеть, значение f (X) которой максимально. Если это значение f (X) больше значения f (X0) исходной (для данной итерации) сети, то из исходной сети убирается узел, благодаря исключению которого и достигается это увеличение значения f (X). Исключение переменных повторяется до тех пор, пока возможно увеличение значения f (X) путем исключения переменных. В качестве значения f (X) нами выбрана величина AUC [8]. Особенно просто проводить оптимизацию БС с наивной топологией и конечной точкой в корневом узле, поскольку удаление любого количества узлов-листьев не меняет топологию, и процесс вычисления AUC точен и быстр.

Алгоритм оптимизации числа узлов БС с наивной топологией реализован разработанной нами программой SiLVIA (Simple Learn Variable Influence Analyzer).

С помощью программы NetGen были построены наивные БС для двух исследуемых К.Т. Построенные сети содержали все выбранные параметры базы данных и обозначались как БС 32+1, где число 32 обозначало количество узлов-листьев. Численная оценка качества прогноза сетей осуществлялась программой ANN. Значение AUC для КТ1-П было равным 0,6263, а для КТ2-С — 0,6209. Эти оценки свидетельствуют о низком качестве прогноза построенных сетей. Далее с помощью программы SiLVIA была осуществлена оптимизация начальных сетей по числу узлов. В связи с тем что количество узлов БС было небольшим (32), программа SiLVIA запускалась с опцией уменьшения на каждом шаге итерации одной переменной до тех пор, пока в сети есть листья. После проведенного анализа результатов оценки качества БС на всех шагах оптимизации для каждой КТ были отобраны по 7 наиболее значимых параметров. Дальнейшее уменьшение числа узлов приводило только к уменьшению AUC. С помощью программы NetGen были построены оптимизированные сети, содержащие только отобранные параметры. На рис. 5 представлены полученные оптимизированные сети.

Результаты оценки качества оптимизированных сетей оказались намного лучше: для КТ1-П было получено значение AUC 0,8331, а для КТ2-С — значение AUC 0,9073. На рис. 1 построены ROC-кривые для различных Б.С. Начальные сети обозначаются как БС 32+1, а оптимизированные — как БС 7+1. На рисунке наглядно видно улучшение прогнозируемой способности оптимизированных БС, особенно сильно увеличилась площадь под ROC-кривой для прогноза смерти пациенток.

Проанализировав оптимизированные сети, мы выявили критические параметры, определяющие риск возникновения КТ1-П и КТ2-С. Таким образом, для прогноза прогрессирования РМЖ (КТ1-П) оказались важны «классические» стандартные факторы риска: гормональный статус опухоли (рецепторы прогестерона), экспрессия рецепторов HER-2/neu и категория N (поражение лимфатических узлов). Известно [27], что гормональный статус опухоли является одним из важнейших прогностических факторов разных исходов РМЖ. Те женщины, в опухоли которых выражена экспрессия гормональных рецепторов, имеют более благоприятный прогноз и меньший риск прогрессирования заболевания. Мы считаем, что совершенно закономерно в число значимых параметров включена экспрессия рецепторов HER-2/neu. Многими работами [28] показано, что опухоли, гиперэкспрессирующие рецепторы HER-2/neu, являются более злокачественными, и пациентки этой группы прогрессируют значительно раньше. Категория N в стадии означает вовлечение регионарных лимфатических узлов в опухолевый процесс. Этот фактор также является значимым для прогноза риска прогрессирования [28].

В нашем исследовании оказалось, что возраст на момент операции вошел в число значимых факторов для КТ1-П. Эти данные согласуются с большим количеством работ, в которых показано, что женщины в возрасте 30—40 лет имеют худший прогноз, чем женщины старше 60 лет, при этом разница в длительности безрецидивного периода составляет 4 года [29]. Биологический возраст, т. е. менструальный статус, является значимым фактором как для прогноза прогрессирования, так и для смерти больных. Известно [29], что у молодых менструирующих женщин злокачественный процесс в молочной железе протекает более агрессивно, и эта когорта пациенток быстро прогрессирует.

Из вариантов лечения значимыми для прогноза прогрессирования РМЖ стали предоперационная лучевая терапия (обсуждается ниже) и неоадъювантная химиотерапия. С одной стороны, назначение неоадъювантной терапии является следствием наличия РМЖ Т3—Т4 стадии, т. е. неблагоприятным прогностическим признаком, с другой стороны — в этой группе неоадъювантная химиотерапия увеличивает безрецидивную выживаемость.

Для прогнозирования смерти больных (КТ2-С) оказались значимы как стандартные прогностические факторы (категория Т, категория N) и молекулярно-биологический подтип опухоли, так и новый прогностический маркер (YB-1). Его повышенная экспрессия ассоциирована с пролиферативным и метастатическим потенциалом опухолевых клеток, и является фактором риска более быстрого наступления смерти больных.

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

Для обеих КТ (прогрессирования и смерти) значимым оказался такой вид лечения, как предоперационная лучевая терапия. Этот вид лечения более не применяется в онкологических клиниках в связи с неэффективностью. Предоперационную лучевую терапию в нашей выборке проводили только 26 (8%) больным, которые получали лечение в 90-е годы прошлого века. Из них 9 (35%) пациенток умерли.

Итак, мы построили и обучили БС с достаточно высоким качеством прогноза. Возникает вопрос: как их использовать для прогноза риска наступления соответствующих событий у конкретной больной? Для этих целей удобно использовать гистограммы риска, связывающие условную вероятность состояния КТ, полученную при опросе обученной БС, с априорной вероятностью прогнозирования для пациенток данной группы риска, т. е. с частотой заданных событий в данной группе риска. Гистограммы риска, представленные на рис. 6, были построены на основании результатов, полученных при построении ROC-кривых (см. рис. 1). Условные вероятности мы разбили на четыре одинаковых интервала, каждый из которых соответствует определенной группе риска. Их можно использовать, чтобы определить вероятность наступления события у новой пациентки. Например, при опросе обученной БС условная вероятность состояния КТ1-П со значением 1 для новой пациентки оказалась в группе риска с диапазоном от 0,25 до 0,5. В этом случае мы можем рассмотреть соотношение известных исходов данного события для пациенток из этой группы риска нашей базы данных. Там оказались 36 пациенток, из которых состояние КТ1-П со значением 1 было у 4. Это означает, что вероятность прогрессирования заболевания у новой пациентки будет равна 11% (4/36).

Построенные гистограммы показывают врачу риск наступления события у пациентов определенных групп. Новый пациент попадает в ту или другую группу риска на основании опроса по его данным обученной Б.С. Например, определим вероятность прогноза прогрессии заболевания, если при опросе нового пациента на обученной БС условная вероятность КТ1-П была равной 0,8. На левой гистограмме изображено 4 группы риска для КТ1-П. По горизонтальной оси отложена условная вероятность наступления события, на основании которой мы определим группу риска для нового пациента, это будет 4-я группа от 0,75 до 1. Далее по вертикальной оси мы определим вероятность прогноза прогрессии заболевания для 4-й группы, которая будет равна 53%. Заметим, что условная вероятность состояния КТ может отличаться от самой вероятности прогноза.

Дополнительно было проведено аналогичное исследование БС для прогноза КТ с момента окончания противоопухолевого лечения, т. е. срок наблюдения сдвигался в среднем на полгода: КТ3-П для прогноза прогрессии заболевания пациента и КТ4-С для прогноза смерти пациента. Временной интервал также составлял 4 года. Результаты оценки качества оптимизированных сетей незначительно различались от первоначального исследования: для КТ3-П значение AUC равнялось 0,8267, а для КТ4-С — 0,9046. Следует отметить, что узлы оптимизированных сетей КТ2-С и КТ4-С получились одинаковыми, а узлы оптимизированных сетей КТ1-П и КТ3-П довольно сильно различаются. В оптимизированной сети прогнозирования КТ3-П по сравнению с оптимизированной сетью прогноза КТ1-П оказался критически важным узел, определяющий размер опухоли, — стадия Т. Размер опухоли является значимым параметром, у пациентов с крупными опухолями прогрессирование и смерть в большинстве случаев наступают в более ранние сроки.

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

Таким образом, показана эффективность использования БС с наивной топологией для выявления критических для прогноза течения РМЖ параметров пациенток, найдены оптимальные БС с высоким качеством прогнозирования, а также построены гистограммы риска, которые можно использовать в экспертной системе персонифицированной медицины.

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

Данная работа выполнена в рамках проекта МГУ им. М.В. Ломоносова «Математическое и молекулярное моделирование для медицины и фармакологии».

Выражаем благодарность профессору А.А. Ставровской и академику Л.П. Овчинникову за научную и клиническую помощь на всех этапах данной работы.

Участие авторов:

Концепция и дизайн исследования: Г. П.Г., В.Б.С.

Сбор и обработка материала: Г. П.Г., О.Г.О., Н.И.М., Е.Ю.Р.

Статистическая обработка данных: А.В.С., И.А.С., В.Б.С.

Написание текста: Г. П.Г., А.В.С., Н.И.М., В.Б.С.

Редактирование: Г. П.Г., О.Г.О., Н.И.М., Л.З.В., И.И.С., В.Б.С.

Конфликт интересов: отсутствует.