Кристаллография, 2021, T. 66, № 5, стр. 730-737

ВВЕДЕНИЕ

Подавляющее большинство систем в природе, представляющих интерес для научных исследований в различных отраслях, в частности медицине, биологии и промышленности, являются многокомпонентными и полидисперсными. Описание и расшифровка структуры такого рода систем, не исследованных ранее, в настоящее время являются актуальной задачей. Одним из эффективных методов по восстановлению структуры исследуемого объекта является метод малоуглового рассеяния (МУР) [1–5]. Для анализа полидисперсных объектов в растворе по данным МУР разработан ряд программ, таких как MIXTURE [1], GNOM [6], SASFIT [7], SASVIEW [8], McSAS [9] и другие.

Ранее для оценки качества работы программы MIXTURE (с реализованным методом переменной метрики в варианте Бройдена–Флетчера–Голдфарба–Шанно с простыми ограничениями на параметры [10], от англ. Broyden–Fletcher–Goldfarb–Shanno algorithm, BFGS) было проведено исследование на устойчивость решения задачи определения структурных параметров наночастиц для теоретических данных МУР, рассчитанных от различных модельных систем [11, 12]. Использовались системы, состоящие из сферических частиц двух типов, имеющих перекрывающиеся и неперекрывающиеся объемные распределения по размерам. Кроме того, к наборам данных добавлялся шум двух типов (гауссовский и пуассоновский) с различным относительным уровнем и исследовалась степень отклонения найденного решения от заданной теоретической модели. Было показано, что диапазон области стартовых значений параметров, позволяющий восстановить истинное решение, зависит от многих факторов, в частности от относительного вклада параметров компонентов смеси в общую кривую рассеяния и от уровня и вида шума, содержащегося в данных.

На следующем этапе полученный результат проверяли на примере реальной экспериментальной системы – раствора кремнезолей, представляющего собой смесь частиц TM-50 и SM [13]. Предварительно с помощью программы GNOM [6] было выявлено наличие двух типов полидисперсных компонентов. Далее в рамках двухкомпонентной модели сферических частиц исследовалась возможность восстановления решений программой MIXTURE. Установлено, что характеры границ устойчивости решений схожи для теоретических и экспериментальных данных, что свидетельствует о стабильности работы алгоритма, реализованного в программе MIXTURE. В то же время ограничения областей стартовых значений параметров, при которых возможно восстановление, указывает на попадание в локальные минимумы в процессе поиска решения.

Вследствие этого возникла идея проверить эффективность работы модифицированной версии программы MIXTURE, в алгоритме которой реализован метод имитации отжига (англ. Simulated annealing, SA) [14], в ряде случаев позволяющий находить глобальный минимум целевой функции. Результаты такого исследования приведены в [15]. В качестве объекта исследования был выбран тот же раствор кремнезоля, что и в [13]. Анализ сходимости решения проводили при одновременном варьировании только одной пары параметров – средних радиусов частиц кремнезоля R₁ и R₂. В результате было установлено, что обе модификации программы (с использованием BFGS- и SA-процедур) дали схожие результаты, но при этом границы для BFGS-процедуры оказались шире.

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

Цель данной работы – систематическое сравнение границ устойчивого приближения для программы с различными минимизационными схемами. Для этого провели исследования на устойчивость упомянутых выше модификаций программы MIXTURE (BFGS, SA, BFGS + SA) для трех пар параметров: (R₁, R₂), (dR₁, R₂) и (R₁, dR₂), где dR₁ и dR₂ – полидисперсности двух типов частиц кремнезоля. Помимо этого исследовали возможности алгоритма с применением симплекс-метода Нелдера–Мида (англ. Nelder–Mead, NM) [17]. Исследования проводили на примере экспериментальных данных МУР от полидисперсного раствора кремнезоля.

КРЕМНЕЗОЛИ

Раствор, выбранный в качестве объекта исследования, получен путем смешения монодисперсных растворов кремнезолей Ludox TM-50^TM и SM^TM. Образцы ТМ-50 (50 мас. % SiO₂ и 0.3 мас. % Na) и SM (30 мас. % SiO₂ и 0.2 мас. % Na) изготовлены компанией Grace Davidson [18]. Несмотря на то что размеры (диаметры) частиц, заявленные производителем (для сухих образцов), составляют 70 Å для SM и 220 Å для ТМ-50, ранее было показано [19], что в растворе они увеличиваются до 100 Å для SM и 270 Å для ТМ-50. Раствор кремнезоля разводили водой в соотношении 1:100, чтобы избежать эффекта межчастичного взаимодействия, при этом центрифугирование не использовали, так как растворы изначально стабилизировались с использованием гидроксида натрия NaOH для предотвращения агрегации.

ЭКСПЕРИМЕНТАЛЬНАЯ УСТАНОВКА

Эксперимент по исследованию раствора кремнезоля выполнен на установке “АМУР-К” [20] (ФНИЦ “Кристаллография и фотоника” РАН, Москва, Россия) и проведен следующим образом: образец кремнезоля помещали в стеклянный капилляр, а затем в вакуумную камеру. Расстояние образец–детектор составляло 700 мм. Время измерений одного образца составляло 10 мин. Экспериментальные данные нормировали на интенсивность падающего пучка, после чего в них вводили поправку на коллимационные искажения.

Установка “АМУР-К” оснащена однокоординатным позиционно-чувствительным детектором ОД3М при фиксированной длине волны излучения λ, равной 0.1542 нм (CuK_α-линия острофокусной трубки, монохроматор из пиролитического графита), и коллимационной системой Кратки. Сечение рентгеновского пучка составляло 0.2 × 8 мм², область углов рассеяния соответствовала диапазону значений модуля вектора рассеяния

В результате эксперимента получена кривая МУР от раствора кремнезоля (рис. 1а, точки). Наилучшее приближение к экспериментальным данным, рассчитанное программой MIXTURE (рис. 1а, сплошная линия), и объемное распределение по радиусам частиц представлены на рис. 1б. Также приведено приближение, при котором решение не попало в глобальный минимум функции (рис. 1а, пунктирная линия). Согласно эксперименту радиусы равны 49 и 140 Å, что согласуется с результатами в [19]. Параметры модели, восстановленные с использованием четырех модификаций программы MIXTURE, представлены в табл. 1. Видно, что результаты восстановления схожи.

МЕТОД РАСЧЕТА ИНТЕНСИВНОСТИ МАЛОУГЛОВОГО РАССЕЯНИЯ ОТ МНОГОКОМПОНЕНТНОЙ ПОЛИДИСПЕРСНОЙ СИСТЕМЫ

Интенсивность рассеяния многокомпонентной системы в случае разбавленной системы можно представить в виде

где суммирование проводится по разным компонентам, ${{{v}}_{k}}$ – относительная объемная доля k-го компонента, I_k(s) – интенсивность k-го компонента.

Для полидисперсной системы частиц интенсивность каждого компонента можно представить в виде

где R – размер частицы, D_k(R) – нормированное объемное распределение частиц по размеру для k-го компонента, V_k(R) – эффективный объем k-го компонента, ∆ρ_k(R) – рассеивающий контраст k-го компонента, i_0k(s,R) – нормированный формфактор k-го компонента.

Для расчета функции объемного распределения k-го компонента использовали распределение Шульца:

где $z = {{\left( {\frac{{{{R}_{0}}}}{{\Delta R}}} \right)}^{2}} - 1$, R₀ – среднее значение размера частицы, ∆R – дисперсия распределения, Г(z) – гамма-функция.

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

где N – число экспериментальных точек, c – шкалирующий коэффициент, I(s_j) и I_calc(s_j) – экспериментальная и рассчитанная интенсивности соответственно, σ(s_j) – экспериментальная ошибка, соответствующая модулю вектора рассеяния s_j.

Алгоритм, использующий формулы (2)–(4), реализован в программе MIXTURE. Такой подход применим к исследованию систем различной природы, в частности он успешно использовался для количественного описания температурного структурного перехода в трехкомпонентных микроэмульсиях [21], анализа ансамбля состояний самособирающихся икосаэдрических капсул фермента люмазинсинтазы в нативном состоянии и при наличии мутаций [22], а также изучения олигомерного состава ассоциатов наночастиц оксида железа, внедренных в капсулы вируса мозаики Бром [23].

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

ИСПОЛЬЗУЕМЫЕ МЕТОДЫ МИНИМИЗАЦИИ

В ходе работы проводили исследование на устойчивость алгоритмов четырех модификаций программы MIXTURE. В каждой из которых использовали различные минимизационные схемы: метод BFGS [10], SA [14], комбинация этих методов (BFGS + SA), а также NM [17].

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

Метод SA является стохастическим методом поиска глобального экстремума. Он является наиболее подходящим для минимизации функций с большим числом переменных, что идеально подходит для многих задач МУР. Алгоритм основан на методе случайных испытаний, запускается с произвольной модели, параметры которой затем изменяются случайным образом в соответствии со схемой Метрополиса [24] до тех пор, пока решение не станет соответствовать экспериментальным данным, например, по критерию χ². В отличие от метода квазиньютоновского градиентного спуска, где на каждом следующем шаге принимается решение с меньшим χ² (5), в методе SA с определенной долей вероятности (определяемой так называемой “температурой отжига”) может приниматься худшее решение, это дает возможность программе выходить из локальных минимумов и продолжать спуск. По мере продолжения поиска значение температуры монотонно снижается, с использованием шкалирующего коэффициента 0.93–0.95 после каждых 1000–10 000 вариаций модели. Однако поиск решения с помощью такого алгоритма занимает длительное время (так как требуется достаточно большое количество температурных итераций).

Комбинация BFGS и SA представляет собой поочередное использование этих алгоритмов минимизации. Схема поиска минимума показана на рис. 2. В целях поиска первого (возможно, локального) минимума нелинейной целевой функции (1) проводят старт программы MIXTURE, использующей градиентный метод минимизации, с некоторых значений параметров двухкомпонентной полидисперсной модели N₁ и N₂ (где N – один из структурных параметров R₁, R₂, dR₁ либо dR₂). Далее с оптимизированных методом BFGS значений проводится запуск модифицированной версии MIXTURE с использованием метода имитации отжига для того, чтобы вывести систему «из равновесия», причем со всего лишь одной итерацией по температуре, поэтому вычисления не занимают длительного времени. После этого с новых полученных значений программа запускается вновь с градиентным методом и так повторяется до тех пор, пока параметры модели не перестанут изменяться. За пороговое значение их изменчивости брали 1–5% от величины, полученной на предыдущей итерации.

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

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

ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ ВОССТАНОВЛЕНИЯ РЕШЕНИЙ

Приближение проводили в рамках модели, состоящей из двух типов сферических частиц с радиусами R₁, R₂, полидисперсностями dR₁, dR₂ и объемными долями V₁, V₂. Исследовали сходимость расчетов к правильному решению при одновременном варьировании стартовых значений следующих пар параметров: (R₁, R₂), (dR₁, R₂), (R₁, dR₂), тогда как остальные параметры фиксировались на истинных значениях, найденных программой MIXTURE (табл. 1) с относительными объемными долями V₁ = 0.25 для меньших частиц кремнезоля (SM) и V₂ = 0.75 для больших частиц кремнезоля (TM-50). В табл. 1 эти значения не приведены, так как алгоритмы всех исследуемых модификаций программы MIXTURE устойчивы к изменению объемных долей во всем диапазоне от 0 до 1.

Диапазоны, в пределах которых выбирали стартовые значения параметров, были следующими: для радиусов R₁ 10.0–96.0 Å и R₂ 100.0–500.0 Å, для полидисперсностей dR₁ и dR₂ 0.5–50.0 Å. Диапазоны были выбраны таким образом, чтобы радиусы меньших и бóльших частиц не коррелировали друг с другом. При этом шаг расчетной сетки по R₁ и R₂ составлял 2 Å, шаг сетки по dR₁ и dR₂ – 0.5 Å. Таким образом, для пары параметров (R₁, R₂) были посчитаны 8844 точки, для (dR₁, R₂) – 20 100 точек, для (R₁, dR₂) – 4400 точек. На основе этих данных построены двумерные карты, представленные на рис. 3. Верхняя панель соответствует результатам, полученным модифицированной версией программы MIXTURE с реализованным методом NM, вторая – программе с методом SA, третья – с BFGS, четвертая (нижняя) панель – комбинации методов BFGS и SA. Первый столбец (слева) – для пары параметров (R₁, R₂), второй – для (dR₁, R₂), третий – для (R₁, dR₂). Область светлого цвета соответствует нахождению истинного решения (в случае, когда значения параметров, полученных в результате моделирования программой MIXTURE, находились в пределах 5% от найденных ранее оптимальных значений). Область темного цвета соответствует всем остальным решениям, существенно отличающимся от истинного решения. Пунктирной линией отмечены истинные (оптимальные) значения параметров, восстановленные программой MIXTURE, указанные в табл. 1.

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

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

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

Анализируя двумерные карты, отметим, что метод NM оказался наименее эффективным по сравнению с другими использовавшимися методами. Несмотря на то что его границы устойчивости шире, чем для метода SA, сравнительно легко попасть в область неустойчивости при незначительном отклонении от истинного значения. Для других трех методов можно говорить об имеющихся тенденциях. Для параметров R₁ и R₂ прослеживаются следующие закономерности: стартовые значения допустимо варьировать в пределах 50, 65 и 80% для R₁ и 15, 65, 100% для R₂ от истинного значения для методов SA, BFGS, BFGS + SA соответственно. Во всех случаях при значении R₂ более 300 Å не удалось добиться восстановления данных. То есть стартовое значение данного параметра может не более чем в 2 раза превышать его истинное значение. При этом решения являются устойчивыми к варьированию параметров dR₁ и dR₂ во всем исследуемом диапазоне.

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

ЗАКЛЮЧЕНИЕ

Исследована устойчивость решений программы MIXTURE для двухкомпонентного раствора кремнезоля с использованием следующих методов минимизации: градиентного метода (BFGS), метода имитации отжига (SA), сочетания этих двух методов (BFGS + SA), а также симплекс-метода Нелдера–Мида (NM). Анализ полученных данных позволяет расположить методы в порядке возрастания их эффективности: NM, SA, BFGS, BFGS + SA. Таким образом, алгоритм поиска с последовательным комбинированием методов минимизации BFGS и SA для систем сферических частиц с малоперекрывающимися объемными распределениями частиц по размерам позволяет максимально эффективно восстанавливать функцию ${{D}_{{v}}}(r)$ в более широких диапазонах стартовых значений параметров, чем с использованием других предложенных алгоритмов. Развитые подходы применимы к исследованию систем различной природы, в том числе для растворов белковых комплексов.

Работа выполнена при поддержке Министерства науки и высшего образования РФ в рамках выполнения работ по Государственному заданию ФНИЦ “Кристаллография и фотоника” РАН, a также при частичной финансовой поддержке Российского фонда фундаментальных исследований (грант № 19-32-90190 Аспиранты) и Программы деятельности НИЦ КИ на 2018–2022 г. Малоугловые измерения выполнены с использованием оборудования ЦКП ФНИЦ “Кристаллография и фотоника” РАН при поддержке Министерства науки и высшего образования России (проект RFMEFI62119X0035).