CC BY
48
УДК 519.71 ] .3:541.12.012.3
П. А. Гуриков, А. В. Колнооченко, Н. В. Меньшутина
Российский химико-технологический университет им. Д.И. Менделеева, Москва. Россия
МОДЕЛИРОВАНИЕ СТРУКТУР АЭРОГЕЛЕЙ И ПРОЦЕССОВ ВЫСВОБОЖДЕНИЯ АКТИВНЫХ ВЕЩЕСТВ
Aerogels' structures model for computer simulations had presented. Cellular automate with Margolœ neiborhood vvas used for describing dynamic ofdrug release. Modelitig algor'ithms were suggested. Kinetic curves of the release of modeling drug were obtained.
В настоящей работе предлагается модель структуры, имитирующей структуру реального аэрогеля. Динамика высвобождения молекул лекарственного (активного) вещества из пористой матрицы аэрогеля в раствор описывается клеточными автоматами с окрестностью Марголуса. Для рассмотренных в работе моделей предложены алгоритмы их реализации. Результатом работы являются модельные кривые высвобождения активного компонента из аэрогеля.
Многие процессы при физическом моделировании описываются дифференциальными уравнениями в частных производных. Методы численного решения таких уравнений давно разработаны и для них существуют эффективные алгоритмы. Однако зачастую простые уравнения имеют сложные граничные условия, что представляет серьезную трудность при их решении. Такие условия часто встречаются при описании массо- и теплопере-носа, динамики биологических популяций, движений толпы в чрезвычайных ситуациях и многих других задачах.
Последние годы ознаменовались разработкой новых методов для описания динамических систем различной природы, которые не основаны на решении дифференциальных уравнений. Такими инструментами являются клеточные автоматы, метод Монте-Карло, теория перколяции и другие.
В настоящей работе клеточными автоматами описывается динамика высвобождения молекул лекарственного (активного) вещества из пористой матрицы в раствор. Это позволяет прогнозировать концентрацию активных веществ в организме и добиваться наиболее значительного терапевтического эффекта. Задачи, аналогичные рассмотренным, возникают при моделировании сушки пористых тел (древесина, керамика), набухании капиллярных тел, процессов адсорбции-десорбции в химической технологии.
Важной особенностью представленной работы является то, что в ней разработан и реализован алгоритм высокопроизводительных параллельных вычислений для метода клеточных автоматов с окрестностью Марголуса с применением новейшей технологии CUDA.
Постановка задачи
В работе [1] в качестве пористой матрицы-носителя изучаются аэрогели - высокопористые тела на основе диоксида кремния с чрезвычайно низкой плотностью (до 0,001 г/см3). Аэрогели являются нетоксичными, процесс их получения сравнительно дешев.
Однако, связь между кинетическими характеристиками высвобождения, структурой аэрогеля и химическим строением активного вещества до сих пор не установлена и носит в лучшем случае эмпирический характер.
В настоящей работе предлагается модель структуры с заданной степенью пористости, имитирующей структуру реального аэрогеля. Для моделирования высвобождения (десорбции) активного вещества в жидкую среду используется метод клеточных автоматов с окрестностью Марголуса. Для обеих моделей разработаны алгоритмы и их программная реализация.
Модель слабо перекрывающихся сфер
Многочисленными экспериментами установлено, что структура аэрогеля представляет собой систему твердых глобул (шаров) диоксида кремния диаметром 2-10 нм. Пространство между ними образует систему пор. Доля свободного пространства, приходящегося на поры, в зависимости от способа получения аэрогеля составляет от 50 до 95 %. В предложенной модели глобулы диоксида кремния моделируются жесткими шарами постоянного радиуса.
Алгоритм построения структуры состоит из следующих шагов:
1. Создание набора шаров N одинакового диаметра й = 2г, перекрывающихся между собой не более определенной величины у/. Первый и последующие шары задаются тройкой случайных координат их центра. Если вновь созданный шар перекрывается с соседним более чем на у«?, то он перемещается вдоль линии, соединяющей центры атомов до тех пор, пока перекрывание не уменьшится до величины у/ (рис. 1). Всякий раз после добавления нового шара вычисляется объем, занимаемый шарами, и пористость -как доля свободного объема. Критерием перехода к шагу 2 является достижение заданной пористости. Если число попыток добавить новый шар превышает наперед заданное критическое значение Ктх, считается, что достигнута гаютнейшая упаковка, и программа переходит к шагу 2;
Рис. 1. Расстояние между сферами уменьшается до с/ш|п = — у2 при добавлении новой сферы
2. Из сгенерированного набора шаров удаляются произвольно выбранные сферы, Поскольку структура реального аэрогеля представляет собой систему связных шаров, то критерием возможности удаления частицы является сохранение перколяционного кластера в рассматриваемом объеме.
Условие связности проверяется с помощью модифицированного алгоритма маркировки кластеров Хошена-Копельмана. Процесс завершается при достижении заданного значения пористости. В случае если ни один шар нельзя удалить без разрушения перколяционного кластера, а заданная пористость еще не достигнута, пользователю предлагается повторить шаг 2.
Поскольку полученная на шаге 2 структура далее будет использована в методе клеточного автомата, ее необходимо перенести на кубическую сетку. Ячейка, в которой находится центр шара и все ее соседи, отстоящие от нее не более чем на г, помечаются как занятые диоксидом кремния (рис. 2).
В экспериментальных исследованиях высвобождения активных веществ из пористой структуры аэрогелей основной величиной, характеризующей количество введенного вещества А, является массовая доля со(А). Для кубической сетки с шагом к легко вычислить, что величина со(А) связана с объемной долей заполнения гр{А) следующим соотношением:
-«АШ_
<р(А)М(А) + <р{8Ю1)р(ЗЮ2)к:> ' ^
где М(А) - молярная масса вешества А, р(8Юз) - плотность кристаллического диоксида кремния, д>(^Юз) - объемная доля диоксида кремния. Таким образом, задавшись определенной массовой долей, по формуле (1) вычисляется доля ячеек, занятых активным веществом. Активное вещество распределяется случайным образом по свободным ячейкам.
Изложенная выше последовательность действий позволяет сгенерировать структуру аэрогеля заданной пористости, содержащую в своих порах одиночные молекулы активного вещества. Далее такая структура исследуется в отношении кинетики высвобождения с помощью метода клеточных автоматов с окрестностью Марголуса.
^— -Ч
— и \ 1
Рис. 2. Перенесение шара на кубическую сетку (сечение)
Клеточные автоматы с окрестностью Марголуса Клеточные автоматы являются абстрактными динамическими системами с дискретным временем. Моделируемая система разбивается на ячейки, каждая из которых в дискретные моменты времени 1 = 1,2,... меняет свое состояние в зависимости от состояния соседних ячеек. Число возможных состояний ячеек конечно. Правила перехода, по которым изменя-
ются состояния автомата, могут быть детерминистическими или вероятностными.
Сложная динамика системы может быть описана при использовании простых правил перехода. Локальные правила могут в явном или неявном виде учитывать законы сохранения (массы, импульса, энергии) или некоторые фундаментальные принципы (возрастание энтропии в изолированной системе, принцип максимального действия). Особым классом клеточных автоматов являются такие, правила которых в явном виде учитывают микроскопическую обратимость физических систем. Одним из простых автоматов с таким поведением является клеточный автомат е окрестностью Маршяуса.
Для наглядности рассмотрим двумерный вариант такого клеточного автомата. Множество клеток разбивается на множество блоков - конечных, однородных частей (рис. За), В простейшем виде блоки имеют размер 2x2. На следующем шаге разбиение меняют так, чтобы возникало некоторое перекрытие блоков на шагах t„.j и t„ (рис. 36). Именно такую схему называют окрестностью Марголуса. Таким образом, разбиения чередуются, чем достигается однородность пространства и времени.
л Б
Рис. 3. Чередование четной (s) и нечегиой (б) решеток на каждой итерации
Для всех блоков задается единое правило перехода: на каждом временном шаге блок может поворачиваться по часовой, стрелке на тс/2, против часовой стрелки я/2 или оставаться без поворота с вероятностями Р«, Рсс* и Р„ соответственно. В простейшем случае вероятности равны 1/3 или P<nv= Рее» =1/2, Рц » 0. В работе [2] показано, что в последнем случае динамика поведения такого автомата удовлетворяет уравнению диффузии с коэффициентом диффузии D = 1,5.
В настоящей работе описанный выше алгоритм распространен на трехмерный случай и состоит из следующих шагов:
1. Трехмерная сетка разбивается на блоки размерностью 2x2x2 двумя способами - четным и нечетным разбиениями;
2. На текущем временном шаге каждый из блоков четного разбиения поворачивается на я/2 по часовой стрелке, против часовой стрелки или сохраняет свое исходное положение с равной вероятностью (выбирается при помощи генератора случайных чисел). Аналогичная операция проделывает-ся с блоками нечетного разбиения. Блоки, в которых хотя бы одна ячейка оказалась занятой структурой диоксида кремния, не поворачиваются. Тем самым моделируются жесткие стенки пор аэрогедя. Ячейки, содержащие
молекулы активного вещества и вышедшие за границу объема структуры, считаются ушедшими в объем растворителя и в дальнейшем не учитываются.
3. Шаги I. 2 повторяются до тех пор. пока в процессе десорбции из структуры не будет удалено заданное количество молекул вещества, или кривая высвобождения перестанет возрастать. Последнее легко проверить статистически.
Рассмотренный алгоритм требует для своей реализации существенных вычислительных затрат. Созданная для реализации алгоритма программа для параллельных вычислений, реализована на высокопроизводительном компьютере, собранном по технологии CUDA. Такая технология позволяет выполнять программы на графическом процессоре (GPU) видеокарты. В настоящей работе использовался компьютер, состоящий из 4 видеокарт (GeForce GTX 280) с объемом памяти t Гб на каждой карте. Программная реализация алгоритмов выполнена на языке программирования С.
Результаты моделирования
На первом этапе для глобул радиуса г = 4 нм и величины перекрывания у/ = 0.4 сгенерировано три структуры аэрогелей с различным значением пористости ш = 60. 80 и 90 %. Па рис, 4 изображен элементы шготнейшей структура (4а) и структуры с частью удаленных шаров (46). Структура на рис 46 является перколирующей, не смотря на кажущееся отсутствие связности шаров в нижней части рисунка. Это противоречие снимается наличием тороидальных граничных условий.
Рис. 4. Исходная плотно упакованная структура (а), перколирующий кластер (б)
Полученные структуры были перенесены на кубическую сетку размером 529x529 х529 с шагом Ь = I нм. Правило перехода, моделирующее диффузию, давалось заданием вероятностей поворотов Р„ = РССЧ¥ = Р„ = 1/3. Параметры в формуле (1) равны: М(А) = 250 г/моль, <р(ЗЮ2) = 2,2 г/мл. На рис. 5 показана динамика высвобождения активного вещества при разных значениях пористости.
Обращает на себя внимание тот факт, что при низкой пористости (60 %) высвобождение происходит быстрее, чем при со = 80 и 90 %, однако абсолютное количество десорбированного вещества мало (не более 1.5 %).
По-видимому, это связано с тем, что повороты ячеек практически не происходят, что отвечает случаю сильной адсорбции. При значениях пористости ш = 80 и 90 % объемная доля высвобожденного вещества оказывается больше: 20 и 55 % соответственно.
Число итерация
Рис. 5, Кинетика высвобождения активного вещества из пористого теля 500x5011x500 им: 1:« - 90 %, 2: ы 80 %, 3: т - fiO %
Модельное время (число итераций) при практическом моделировании необходимо связать с физическим временем, что достигается сравнением экспериментальной кривой с рассчитанной для какого-либо значения пористости.
Заключение
Рассмотренные выше модели позволяют, используя простые правила перехода, описать диффузию в пористых телах. Предложенный алгоритм клеточного автомата с окрестностью Марголуса позволяет проводить параллельные вычисления, что особенно актуально при решении трехмерных задач массо- и тештопереноса и найдет приложение в многочисленных областях физического и химического моделирования. За рамками обсуждения остались вопросы деградации структуры аэрогеля в растворе {т.н. коллапс аэрогеля) и связи между свойствами активного вещества и кинетикой его высвобождения. Последний вопрос сравнительно легко решается в рамках структурно-топологического описания сольватации и станет предметом дальнейших публикаций.
Библиографические ссылки
1. !. Smirnova [ets.J // Pharmaceutical development and technology. V. 9 (2004), N. 4. PP. 443-452.
2. Мдлииецкий Г. Г., Степанцов М. Е. // Ж. выч. мат. и матем. физ„ 1998. Т. 38. N6. С. 1017 - 1020.
