Кристаллография, 2019, T. 64, № 2, стр. 321-326

Моделирование роста биокристаллов с помощью температурного поля

В. И. Стрелов 1*, В. П. Гинкин 2, И. Ж. Безбах 1

1 Институт кристаллографии им. А.В. Шубникова ФНИЦ “Кристаллография и фотоника” РАН
Москва, Россия

2 Государственный научный центр РФ “Физико-энергетический институт им. А.И. Лейпунского”
Обнинск, Россия

* E-mail: kmikran@spark-mail.ru

Поступила в редакцию 06.03.2018
После доработки 27.03.2018
Принята к публикации 03.04.2018

Полный текст (PDF)

Аннотация

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

ВВЕДЕНИЕ

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

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

В эксперименте [6] с вынужденным конвективным потоком при выращивании тетрагональных кристаллов лизоцима обнаружено сильное уменьшение скорости роста со временем. Авторы [7] выращивали тетрагональные кристаллы лизоцима, подвешенные в восходящем потоке раствора в термосифонной петле. При высоких пересыщениях рост кристаллов прекращался, когда их размер достигал ~0.1 мм. Однако затравочные кристаллы в растворах без вынужденных потоков продолжали расти до больших размеров. Эти явления говорят о том, что механизм влияния конвекции на рост кристаллов изучен недостаточно хорошо.

В [3] изучался механизм ступенчатого роста от края грани, а время зарождения грани представлено в вероятностной форме. Ступенчатый рост грани наблюдается на фотографиях в эксперименте и его можно сравнить с расчетным процессом. В [8] для моделирования зарождения и ступенчатого роста биокристаллов использовался метод Монте-Карло. В экспериментальных работах [9, 10] исследовались состояния молекул лизоцима на стадии начала процесса кристаллизации и показано, что в растворе белка присутствуют молекулы олигомеров, кроме того, тип растворителя влияет на начальную стадию кристаллизации, что, конечно, требует уточнения моделей процесса зарождения и роста кристаллов белков.

В [1113] рассматривали объемное осаждение белка из пересыщенного раствора при диффузии осаждающей соли в камеру с белком из прилегающего объема с раствором соли. Обычная геометрия – капилляр, частично заполненный белковым раствором (возможно, белок распределен в геле для подавления конвекции), который находится в контакте с частью, заполненной раствором соли (NaCl). Наблюдали за числом и размером получающихся кристаллов в зависимости от концентраций и пересыщений, меняющихся по длине капилляра. Решали безразмерные уравнения сохранения массы, количества движения и диффузии концентраций белка и соли.

Анализ литературных данных позволяет сделать вывод, что температуру можно использовать как средство управления кристаллизацией белков. В [14] исследованы зависимости растворимости S от температуры для двух типов белков на примерах лизоцима и альбумина лошадиной сыворотки. На основе этих данных в [15] получена двухпараметрическая зависимость растворимости лизоцима S от концентрации осадителя и температуры (рис. 1).

Рис. 1.

Двумерное распределение растворимости S лизоцима в зависимости от температуры T и концентрации осадителя CNaCl.

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

ОПИСАНИЕ ЭКСПЕРИМЕНТА

Тонкостенный капилляр внутренним диаметром 0.7 мм с раствором лизоцима помещается в термостат с постоянной в течение всего времени эксперимента температурой Т0 = 32°С. Начальная концентрация лизоцима в растворе C1(0) = = 50 мг/мл, осадителя NaCl C2(0) = 4 мас. %, т.е. C2(0) = 0.04ρ, где ρ – плотность раствора. В начальный момент времени в центральной части капилляра на участке стенки длиной 0.05 см путем подвода к ней капилляра тонкой иглы, охлаждаемой до температуры 6°C, на стенке капилляра создается область с пониженной температурой Т1 = 6°С. В результате того, что растворимость лизоцима в области раствора вблизи контакта холодной иглы со стенкой капилляра понижается, в ней возникает зародыш кристалла.

В этот момент, чтобы уменьшить скорость роста образовавшегося кристалла, температуру проводника тепла увеличивают до Т2 = 12°C со скоростью, доступной для аппаратуры, и далее не меняют в течение всего эксперимента (приблизительно 6 сут).

На рис. 2а показан зародыш кристалла лизоцима через одни сутки после начала эксперимента, на рис. 2б выросший кристалл через 6 сут в конце эксперимента. Этот эксперимент продемонстрировал принципиальную возможность управляемого выращивания единичных биокристаллов большого размера и высокого структурного качества. Отметим, что момент времени и место появления первичного зародыша кристалла в области раствора с пересыщением, значительно превышающем единицу, точно непредсказуемы (в каком месте в области иглы), так как являются по своей природе вероятностными, а не детерминированными величинами. Так, в эксперименте при температуре Т1 = 6°С примерно через одни сутки был визуально зафиксирован зародыш кристалла немного сбоку от места расположения острия иглы, в других случаях зародыши образовывались точно на острие иглы.

Рис. 2.

Фотографии зародыша кристалла лизоцима в капилляре (внутренний диаметр капилляра 0.7 мм) в начале эксперимента (а) и выросшего кристалла лизоцима в конце эксперимента (б).

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

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

Зарождение кристалла моделируется детерминистически. В каждый момент времени рассчитывается распределение концентрации лизоцима в растворе C1 и вычисляется пересыщение $\sigma = \frac{{{{C}_{1}}}}{S}$, где S – растворимость.

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

В области кристалл–раствор решаются уравнения конвективного тепломассопереноса в приближении Буссинеска:

(1)
$\begin{gathered} \frac{{\partial {\mathbf{V}}}}{{\partial t}} + ({\mathbf{V}}\nabla ){\mathbf{V}} = - \frac{1}{\rho }\nabla p + \nu \Delta {\mathbf{V}} - \\ - \;{{\beta }_{C}}_{{\text{1}}}{\mathbf{g}}({{C}_{{\text{1}}}} - {{C}_{{{\text{1(0)}}}}}) - {{\beta }_{C}}_{{\text{2}}}{\mathbf{g}}({{C}_{{\text{2}}}} - {{C}_{{{\text{2(0)}}}}}) - \\ - \;{{\beta }_{T}}{\mathbf{g}}(T - {{T}_{{\text{0}}}}), \\ \end{gathered} $
(2)
$\operatorname{div} {\mathbf{V}} = 0,$
(3)
$\frac{{\partial T}}{{\partial t}} + {\mathbf{V}}\nabla T = \frac{1}{{{{c}_{p}}\rho }}\nabla \lambda \nabla T,$
(4)
$\frac{{\partial {{C}_{1}}}}{{\partial t}} + {\mathbf{V}}\nabla {{C}_{1}} = \nabla {{D}_{1}}\nabla {{C}_{1}} + {{S}_{1}},$
(5)
$\frac{{\partial {{C}_{2}}}}{{\partial t}} + {\mathbf{V}}\nabla {{C}_{2}} = \nabla {{D}_{2}}\nabla {{C}_{2}} + {{S}_{2}},$
где ${{C}_{1}}$, ${{C}_{2}}$ – концентрации лизоцима и осадителя (NaCl), ${{C}_{{1(0)}}}$, ${{C}_{{2(0)}}}$ – начальные концентрации лизоцима и осадителя, ${{D}_{1}}$, ${{D}_{2}}$ – коэффициенты диффузии лизоцима и осадителя, V вектор скорости течения раствора, р – давление, h – удельная энтальпия, T – температура, ${{T}_{0}}$ – начальное значение температуры, ${{c}_{p}}$ – теплоемкость, βТ – коэффициент теплового расширения раствора, β1, β2 – коэффициенты концентрационного расширения белка и осадителя, g вектор остаточной гравитации, $\rho = {{\rho }_{{{{{\text{H}}}_{{\text{2}}}}{\text{O}}}}}(1 + {{\beta }_{1}}{{C}_{1}} + {{\beta }_{2}}{{C}_{2}})$ – плотность раствора, ${{\rho }_{L}}$ – плотность кристалла, ν – кинематическая вязкость, ${{S}_{1}}$, ${{S}_{2}}$ – стоки концентраций примесей в кристалл при его росте (либо источники концентраций примесей при его растворении).

Выражения для источников ${{S}_{1}}$, ${{S}_{2}}$ в уравнениях (4), (5) взяты из [16] по аналогии с моделью кристаллизации расплавов в методе направленной кристаллизации Бриджмена и имеют вид:

(6)
${{S}_{1}} = \frac{{\partial {{\varepsilon }_{S}}}}{{\partial t}}(1 - {{K}_{{p1}}}){{C}_{1}} + {{\varepsilon }_{S}}\frac{{\partial {{C}_{1}}}}{{\partial t}},$
(7)
${{S}_{2}} = \frac{{\partial {{\varepsilon }_{S}}}}{{\partial t}}(1 - {{K}_{{p1}}}){{C}_{2}} + {{\varepsilon }_{S}}\frac{{\partial {{C}_{2}}}}{{\partial t}},$
где ${{\varepsilon }_{S}} = \frac{{{{V}_{S}}}}{{{{V}_{{Cell}}}}}$ и ${{V}_{S}}$ – доля и объем твердой фазы в ячейке сетки, ${{V}_{{Cell}}} = {{V}_{S}} + {{V}_{L}}$ – объем ячейки сетки, ${{V}_{L}}$ – объем жидкой фазы в ячейке сетки, ${{K}_{{p1}}}$ и ${{K}_{{p2}}}$ – равновесные коэффициенты распределения белка и осадителя.

Производная от доли твердой фазы по времени $\frac{{\partial {{\varepsilon }_{S}}}}{{\partial t}}$ в выражениях (6), (7) означает скорость движения границы фазового перехода. Если она положительна, то кристалл растет, если отрицательна – растворяется.

Подставив (6), (7) в (4), (5) и проведя несложные преобразования, получим вместо (4), (5) следующие уравнения:

(8)
${{\varepsilon }_{L}}\frac{{\partial {{C}_{1}}}}{{\partial t}} + {\mathbf{V}}\nabla {{C}_{1}} = \nabla {{D}_{1}}\nabla {{C}_{1}} + \frac{{\partial {{\varepsilon }_{S}}}}{{\partial t}}(1 - {{K}_{{p1}}}){{C}_{1}},$
(9)
${{\varepsilon }_{L}}\frac{{\partial {{C}_{2}}}}{{\partial t}} + {\mathbf{V}}\nabla {{C}_{2}} = \nabla {{D}_{2}}\nabla {{C}_{2}} + \frac{{\partial {{\varepsilon }_{S}}}}{{\partial t}}(1 - {{K}_{{p2}}}){{C}_{2}}.$
На фиксированной сетке значения ${{C}_{1}}$ и ${{C}_{2}}$ в узлах сетки являются усредненными по объему ячейки сетки концентрациями лизоцима и осадителя соответственно. Коэффициент ${{K}_{p}}$ равен отношению истинных концентраций примеси в твердой и жидкой фазах на границе фазового перехода. Для осадителя NaCl ${{K}_{{p2}}}$ = 0.01. Для лизоцима ${{K}_{{p1}}} = \frac{{{{\rho }_{L}}{{\varepsilon }_{L}}}}{{{{C}_{1}}}}$. Подставив это выражение в уравнение (8) и проведя несложные преобразования, получим следующее уравнение:
(10)
$\frac{{\partial ({{\varepsilon }_{L}}{{C}_{1}})}}{{\partial t}} + {\mathbf{V}}\nabla {{C}_{1}} = \nabla {{D}_{1}}\nabla {{C}_{1}} - \frac{{\partial {{\varepsilon }_{S}}}}{{\partial t}}{{\rho }_{L}}{{\varepsilon }_{L}}.$
Полученное уравнение эквивалентно (8), но удобней для решения, так как не содержит неизвестное заранее значение ${{K}_{{p1}}}$.

Система пяти уравнений (1)–(3), (9), (10) содержит шесть неизвестных ${\mathbf{V}},\;p,\;T,\;{{C}_{1}},\;{{C}_{1}},\;{{\varepsilon }_{L}}$. Для замыкания уравнений (1)–(3), (9), (10) необходимо сформулировать условие для вычисления ${{\varepsilon }_{L}}$. Это условие получается из требования равенства концентрации лизоцима ${{C}_{1}}$ растворимости $S$ в ячейках сетки с $0 < {{\varepsilon }_{L}} < 1$ и имеет вид

(11)
$\frac{{{{C}_{1}}}}{{{{\varepsilon }_{L}}}} = S.$
Система уравнений (1)–(3), (9)–(11) дополняется граничными и начальными условиями и решается неявным методом. Отметим, что в ячейках сетки с $0 < {{\varepsilon }_{L}} < 1$ все теплофизические данные, входящие в уравнения (1)(3), (9)(11), вычисляются с массами жидкой и твердой фаз.

В этой модели емкость с раствором рассматривается как замкнутая прямоугольная область высотой 1 мм и шириной 6 мм. Для температуры на левой границе области при x = 0 задавалось условие симметрии: $\frac{{\partial T}}{{\partial x}} = 0$. На остальных границах области температура задавалась равной $T = {{T}_{0}}$ = = 32°C, за исключением небольшого участка по оси x вблизи начала координат, где задавалась пониженная температура, а именно: при 0 < x < a, y = 0, $T = {{T}_{2}} = 12^\circ {\text{C}}$ (а = 0.025 см).

Для концентраций ${{C}_{1}}$ и ${{C}_{2}}$ на всех границах жидкой и твердой фаз и на оси симметрии задавалось условие $\frac{{\partial {{C}_{i}}}}{{\partial n}} = 0$, $(i = 1,2)$, где n – нормаль к границе. На правой границе области также задавалось условие $\frac{{\partial {{C}_{i}}}}{{\partial x}} = 0$.

Сетка задавалась равномерной. По оси x число узлов было равно 500, по оси y – 50. Скорость конвекции полагалась равной нулю, что соответствует нулевой силе тяжести. Теплофизические данные для расчета взяты из [15].

АЛГОРИТМ РЕШЕНИЯ

Вводится дискретизация по временной координате с шагом ${{\tau }_{i}}$, где i – номер шага по времени. Уравнения (1)(3), (9)(11) переходят в следующие уравнения:

(12)
$\begin{gathered} \frac{{{{{\mathbf{V}}}^{i}} - {{{\mathbf{V}}}^{{i - 1}}}}}{{{{\tau }_{i}}}} + ({{{\mathbf{V}}}^{i}}\nabla ){{{\mathbf{V}}}^{i}} = - \frac{1}{{{{\rho }^{i}}}}\nabla {{p}^{i}} + \nu \Delta {{{\mathbf{V}}}^{i}} - \\ - \;{{\beta }_{C}}_{{\text{1}}}{\mathbf{g}}(C_{1}^{i} - {{C}_{{{\text{1(0)}}}}}) - {{\beta }_{C}}_{{\text{2}}}{\mathbf{g}}(C_{2}^{i} - {{C}_{{{\text{2(0)}}}}}) - \\ - \;{{\beta }_{T}}{\mathbf{g}}({{T}^{i}} - {{T}_{{\text{0}}}}), \\ \end{gathered} $
(13)
$\operatorname{div} {{{\mathbf{V}}}^{i}} = 0,$
(14)
$\frac{{{{T}^{i}} - {{T}^{{i - 1}}}}}{{{{\tau }_{i}}}} + {{{\mathbf{V}}}^{i}}\nabla {{T}^{i}} = \frac{1}{{c_{p}^{i}{{\rho }^{i}}}}\nabla {{\lambda }^{i}}\nabla {{T}^{i}},$
(15)
$\frac{{\varepsilon _{L}^{i}C_{1}^{i} - \varepsilon _{L}^{{i - 1}}C_{1}^{{i - 1}}}}{{{{\tau }_{i}}}} + {{{\mathbf{V}}}^{i}}\nabla C_{1}^{i} = \nabla D_{1}^{i}\nabla C_{1}^{i} - \frac{{\partial \varepsilon _{S}^{i}}}{{\partial t}}\rho _{L}^{i}\varepsilon _{L}^{i},$
(16)
$\varepsilon _{L}^{i}\frac{{\partial {{C}_{2}}^{i}}}{{\partial t}} + {{{\mathbf{V}}}^{i}}\nabla C_{2}^{i} = \nabla D_{2}^{i}\nabla C_{2}^{i} + \frac{{\partial \varepsilon _{S}^{i}}}{{\partial t}}\left( {1 - {{K}_{{p2}}}} \right)C_{2}^{i},$
(17)
$\frac{{C_{1}^{i}}}{{\varepsilon _{L}^{i}}} = {{S}^{i}}.$

Система уравнений (12)–(17) на каждом шаге по времени решается неявным методом установления. Этот метод заключается в следующем.

Вводится вспомогательная сетка по времени установления с шагом $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\tau } $. К левым частям уравнений (12), (14)–(16) прибавляются разностные производные по времени на этой временной сетке от искомых функций, уравнение для давления получается из разностного уравнения неразрывности, а уравнение (17) заменяется на следующее:

(18)
$\varepsilon _{L}^{i} = \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } {{_{L}^{i}}_{L}} - \frac{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } {{{_{L}^{i}}}_{L}}}}{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\rho } {{{_{L}^{i}}}_{L}}}}\left( {\frac{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{C} _{1}^{i}}}{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } {{{_{L}^{i}}}_{L}}}} - {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S} }}^{i}}} \right).$

Полученная система уравнений на каждом шаге по времени решается до установления по следующему алгоритму.

Шаг 1. По заданному распределению ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }_{L}}$ на предыдущем шаге установления вычисляются ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{C} }_{1}}$, ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{C} }_{2}}$, $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{T} $, $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S} $ и распределение пересыщений $\sigma $ по формуле

(19)
$\sigma = \left\{ \begin{gathered} 0,\quad {\text{п р и }}\quad {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }}_{L}} = 0, \hfill \\ \frac{{{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{C} }}_{1}}}}{{{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }}_{L}}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S} }},\quad {\text{п р и }}\quad 0 < {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }}_{L}} < 1. \hfill \\ \end{gathered} \right.$

Шаг 2. Рассчитывается распределение ${{\tilde {\varepsilon }}_{L}}$ по следующему алгоритму:

если (${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }_{L}} = 1$ и $\sigma < \eta $), то ${{\tilde {\varepsilon }}_{L}} = {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }_{L}}$,

если (${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }_{L}} = 1$ и $\sigma \geqslant \eta $) или $0 < {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }_{L}} < 1$, то

(20)
${{\tilde {\varepsilon }}_{L}} = {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }_{L}} - \frac{{{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }}_{L}}}}{{{{\rho }_{L}}}}\left( {\frac{{{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{C} }}_{1}}}}{{{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }}_{L}}}} - \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S} } \right).$

Шаг 3. Рассчитывается распределение ${{\varepsilon }_{L}}$ по следующему алгоритму.

Если в какой-нибудь ячейке получилось ${{\tilde {\varepsilon }}_{L}} < 0$, то это означает, что вся ячейка заполнилась кристаллом и еще осталась доля твердой фазы $\delta = \left| {{{{\tilde {\varepsilon }}}_{L}}} \right|$, которая должна перейти в одну из соседних ячеек. Положим, что эта избыточная доля твердой фазы переходит в ту ячейку, в которой доля жидкой фазы ${{\tilde {\varepsilon }}_{L}}$ имеет наибольшее значение. Тогда в данной ячейке полагаем ${{\varepsilon }_{L}} = 0$, а в соседней ячейке с наибольшей долей жидкой фазы ${{\tilde {\varepsilon }}_{L}}$ полагаем ${{\varepsilon }_{L}} = {{\tilde {\varepsilon }}_{L}} - \delta $. Если в какой-нибудь ячейке получилось ${{\tilde {\varepsilon }}_{L}} > 1$, то это означает, что весь кристалл в данной ячейке растворился и еще осталась доля жидкой фазы $\delta = {{\tilde {\varepsilon }}_{L}} - 1$, которая должна перейти в одну из соседних ячеек. Положим, что эта избыточная доля жидкой фазы переходит в ту ячейку, в которой доля жидкой фазы ${{\tilde {\varepsilon }}_{L}}$ имеет наименьшее значение. Тогда в данной ячейке полагаем ${{\varepsilon }_{L}} = 1$, а в соседней ячейке с наименьшей долей жидкой фазы ${{\tilde {\varepsilon }}_{L}}$ полагаем ${{\varepsilon }_{L}} = {{\tilde {\varepsilon }}_{L}} + \delta $.

В остальных ячейках полагаем ${{\varepsilon }_{L}} = {{\tilde {\varepsilon }}_{L}}$.

Шаг 4. Сравниваются ${{\varepsilon }_{L}}$ и ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }_{L}}$. Если нужная точность не достигнута, то полагаем ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon } }_{L}} = {{\varepsilon }_{L}}$ и возвращаемся на Шаг 1.

РЕЗУЛЬТАТЫ РАСЧЕТОВ

Расчеты проводились при условии отсутствия силы тяжести g = 0. Они показали, что распределение концентрации осадителя ${{C}_{2}}$ быстро устанавливается равным константе, равной ${{C}_{{2(0)}}}$ во всей области раствора. Это происходит потому, что коэффициент диффузии осадителя на порядок больше коэффициента диффузии лизоцима. Вытесненный на первой стадии роста из области выросшего кристалла осадитель NaCl быстро диффундирует в раствор и равномерно распределяется в нем. Поэтому уравнение (5) можно не решать после первых 100 с кристаллизации, полагая ${{C}_{2}} = {{C}_{{2(0)}}}$ в области раствора.

Температурное поле устанавливается очень быстро, за время менее 100 с. Рост кристалла определяется градиентом концентраций, возникшим после зарождения кристалла в начале системы координат. При этом реализуется ступенчатый рост кристалла, наблюдаемый экспериментально. На рис. 3 представлена динамика роста кристалла в разные моменты времени. Из него следует, что к моменту времени t = 600 000 с рост кристалла практически прекращается. Сравнение с экспериментом показало, что в целом процесс зарождения и роста кристаллов по данной модели воспроизводится верно [17].

Рис. 3.

Динамика роста кристалла в различные моменты времени: 20 (а), 200 000 (б), 1 000 000 с после начала процесса кристаллизации (в).

На рис. 4а представлены расчетные зависимости массы кристалла от времени при двух значениях конечной температуры охлаждающей иглы: 12 и 20°C. Из рисунка видно, что чем меньше температура теплоотвода капилляра, тем выше скорость роста и больше масса выросшего кристалла. Однако из опыта известно, что чем выше скорость роста кристалла, тем хуже его качество. Именно по этой причине трудно вырастить кристалл большого размера и высокого качества.

Рис. 4.

Зависимости массы кристалла m от времени t: при двух значениях температуры охлаждающей иглы (а); с учетом (1) и без учета (2) конвекции (б).

Был также выполнен расчет при нормальной силе тяжести, соответствующей земным условиям ($g = {{g}_{0}}$). Рисунок 4б демонстрирует рост массы кристалла с учетом (кривая 1) и без учета конвекции, вызываемой силой тяжести (кривая 2). Из него видно, что на начальном этапе учет конвекции приводит к значительному ускорению роста кристалла, что соответствует экспериментальным наблюдениям [18].

ЗАКЛЮЧЕНИЕ

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

Работа выполнена при поддержке Министерства науки и высшего образования в рамках выполнения работ по Государственному заданию ФНИЦ “Кристаллография и фотоника” РАН.

Список литературы

  1. Lin H., Rosenberger F., Alexander J.I.D. et al. // J. Cryst. Growth. 1995. V. 151. P. 153.

  2. Monaco L.A., Rosenberger F. // J. Cryst. Growth. 1993. V. 129. P. 465.

  3. Vekilov P.G., Lin H., Rosenberger F. // Phys. Rev. E. 1997. V. 55. P. 3202.

  4. Vekilov P.G., Ataka M., Katsura T. // J. Cryst. Growth. 1993. V. 130. P. 317.

  5. Vekilov P.G. // Progr. Cryst. Growth Charact. 1993. V. 26. P. 25.

  6. Pusey M., Witherow W., Naumann R. // J. Cryst. Growth. 1988. V. 90. P. 105.

  7. Nyce T.A., Rosenberger F. // J. Cryst. Growth. 1991. V. 110. P. 52.

  8. Durbin S.D., Feher G. // J. Cryst. Growth. 1991. V. 110. P. 41.

  9. Марченкова М.А., Волков В.В., Благов А.Е. и др. // Кристаллография. 2016. Т. 61. № 1. С. 10.

  10. Бойкова А.С., Дьякова Ю.А., Ильина К.Б. и др. // Кристаллография. 2017. Т. 62. № 6. С. 876.

  11. Lappa M., Castagnolo D., Carotenuto L. // Physica A. 2002. V. 314. P. 623.

  12. Lappa M., Piccolo C., Carotenuto L. // J. Cryst. Growth. 2003. V. 254 P. 469.

  13. Lappa M., Carotenuto L. // Int. J. Microgravity Res. Appl. 2002. V. XIV. P. 2.

  14. Rosenberger F., Howard S.B., Sowers J.W. et al. // J. Cryst. Growth. 1993. V. 129. P. 1.

  15. Гинкин В.П., Ганина С.М., Стрелов В.И. и др. // Поверхность. Рентген., синхротр. и нейтр. исслед. 2009. № 2. С. 17.

  16. Timchenko V., Chen P.Y.P., G. de Vahl Davis et al. // Int. J. Heat Mass Transfer. 2000. № 43. P. 963.

  17. Безбах И.Ж., Захаров Б.Г., Стрелов В.И., Гинкин В.П. // Тр. 6-й междунар. конф. “Рост кристаллов и тепломассоперенос” (ICSC-2005), Обнинск, 25–30 сентября 2005, С. 545.

  18. Стрелов В.И., Захаров Б.Г., Безбах И.Ж. и др. // Кристаллография. 2008. Т. 53. № 1. С. 164.

Дополнительные материалы отсутствуют.