Исследование механизма трансформации гауссовского распределения концентраций при аномальных режимах изоэлектрического фокусирования Текст научной статьи по специальности «Математика»

УДК 004.942

ИССЛЕДОВАНИЕ МЕХАНИЗМА ТРАНСФОРМАЦИИ ГАУССОВСКОГО РАСПРЕДЕЛЕНИЯ КОНЦЕНТРАЦИЙ ПРИ АНОМАЛЬНЫХ РЕЖИМАХ ИЗОЭЛЕКТРИЧЕСКОГО ФОКУСИРОВАНИЯ

© 2012 г Л.В. Сахарова

Ростовский филиал Морской государственной академии, Rostov Branch of Maritime State Academy,

ул. Седова, 8, корп. 3, г. Ростов н/Д, 344006 Sedov St., 8, build. 3, Rostov-on-Don, 344006

Исследована математическая модель создания естественного рН-градиента в водном растворе амфолитов. Задача возникает в практике использования метода изоэлектрического фокусирования (ИЭФ), применяемого для электрофоретического разделения биохимических смесей. Показано, что в случае сверхвысоких плотностей электрического тока взамен классического распределения концентраций возникает «аномальное» распределение: стандартное гауссовское распределение трансформируется в «платообразное». Дана физическая интерпретация полученному результату, сделаны практические выводы о способах повышения разрешающей способности метода ИЭФ.

Ключевые слова: математическая модель, изоэлектрическое фокусирование, амфолиты, Гауссовское распределение.

The mathematical model of creation the natural pH-gradient in the aqueous solution of ampholytes is investigated. The problem is arising in the practice of isoelectric focusing (IEF); it is laboratory method for separating of biochemical mixtures. In case of overcurrent the classical Gaussian distribution is replaced by anomalous plateau-shaped distribution. The physical interpretation of results is given. The practical conclusions about methods of increasing the IEF-resolution are presented.

Keywords: mathematical model, isoelectric focusing, ampholytes, Gaussian distribution.

Изоэлектрическое фокусирование (ИЭФ) относится к одним из наиболее эффективных и универсальных современных методов фракционирования и анализа белков, является неотъемлемой частью экспериментальных исследований практически во всех отраслях современной биологии: биохимии и биотехнологии, в популяционной и молекулярной генетике, при расшифровке первичной структуры генов, при тонком фракционировании белков, близких по своим физико-химическим свойствам.

Важная проблема ИЭФ - повышение разрешающей способности метода. Однако на практике в ряде случаев задача полного разделения амфолитов представляет собой существенную трудность [1] вследствие неустойчивости градиента рН и перепадов электропроводности. Невозможность точного теоретического прогнозирования картины ИЭФ в этих случаях связана с недостаточно четким пониманием механизма процессов, происходящих в электрофоретической колонке (ЭК). В настоящее время серьезные вопросы вызывает теория гауссовского распределения концентраций, изначально заложенная в основе математической модели ИЭФ [1]. Обычно полагают, что распределение концентрации компонент смеси имеет гауссов вид: С = С ехР^_ рЕх2 / 2В^, где С - концентрация; Е - напряженность поля; В - коэффициент диффузии. Практика подтверждает, что гауссово распределение применимо к широкому классу так называемых амфолитов-носителей, т.е. амфотерных веществ, обладающих хорошими проводимостью и буферной емкостью. В недавно опубликованных работах [2, 3] подвергается сомнению универсальность

гипотезы о гауссовском распределении концентраций: представленные профили концентраций имеют деформированный, «прямоугольный» вид.

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

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

Постановка задачи

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

ности / , константы диссоциации реакций K к 1

(k )

к = 1,2,...N.

(к)

К , а также общие количества т ,

2 к'

Температура Т внутри ЭК считается постоянной. Под действием постоянного тока плотности J в ЭК формируется стационарное распределение концентраций амфолитов, которое приводит к стационарному рас-

и

пределению концентрации ионов водорода - образуется так называемый естественный градиент pH .

Предполагается, что в водном растворе диссоциация k -го амфолита протекает по схеме:

NH+RCOOH о NH RCOOH + H+,

3 2

НИ2RCOOH о никоею- + И+.

21

Здесь НИ ЯСООИ, НИ ЯСОО , НИ RCOOH -

3 2 2

положительный, отрицательный и «нейтральный» ионы амфолита, молярные концентрации которых

обозначим С, С , С. Общая, или так называемая 1 -1 о

аналитическая, концентрация амфолита определяется формулой £ = + ^ + £к . Помимо указанных реакций диссоциации амфолитов в водном растворе следует учитывать реакцию автопротолиза воды

И2O о OH -+ И+.

В равновесном состоянии концентрации ионов связаны с аналитической концентрацией Ск, моль/л, следующими формулами:

Кк к К гk к К ¡,к к к с = а с , с = а с , с = (1 — а — а )с , Ь1 1 ь к -1 -Г к ь0 У 1 —Гь к

к _ а1

H

2

К v") Кх") + K(k) H + H 2

1

2

1

к

а =

2 (к )тлк)

К(к) К(к) 1 2

К V) К(к) + К(к) H + H 2 1 2 1

где ак и а^ - степени диссоциации амфолита; И -

концентрация ионов водорода.

Стационарное распределение концентраций описывается следующей системой уравнений:

СС к к

— е--^ + С (а —а )Е = 0, к = 1,2,...Н ; (1)

Ск к 1 —1

N

J = S

к=1

п d (( к к V) ( к к |е „

-D--1 \ а -а \Е \ + и\ а + а \Е E

к dx II 1 -1) к ) к \ 1 2 ) к

- D ■ — + м HE +

H dx H

+ D

d (OH)

OH

dx

+ u ■ OH ■ E ;

OH

N к к S (а -а + H - OH = 0; к=1

1

-1 к

(2)

(3)

2

l

ж (х)Сх = тк, к = 1,2,...Н . (4)

Здесь Е - напряженность электрического поля; г

к 2

и I - радиус ЭК и ее длина; OH = —к— концентрация

И

2 —14

гидроксил-ионов, к = 10 - так называемое ионное

w

произведение воды; / , / - подвижности ионов

И OH

водорода и гидроксил-ионов; Б , Б , Б - коэф-

к И OH

фициенты диффузии ионов; Б =е/ ; е= ЯТ / Г -

к к

стандартный электрохимический параметр, где величины Я , Т и Г - соответственно универсальная газовая постоянная, температура и число Фарадея.

Соотношение (1) - уравнение массопереноса, полученное на основании уравнения потока амфолита и

С

основного уравнения теории переноса —— + = 0

к

(в стационарном случае и при отсутствии потока вещества через границы ЭК); (2) представляет собой обобщенный, т.е. с учетом диффузии, закон Ома (плотность тока является суммой плотностей токов всех ионов, включая ионы водорода и гидроксила); (3) - уравнение электронейтральности, интегральное условие (4) соответствует закону сохранения массы вещества, в соответствии с которым суммарное количество всех

трех форм амфолита неизменно и равно т . Решение

к

задачи (1)-(4) позволяет найти распределение С (х) и

к

рИ (х) .

Основная математическая трудность численного интегрирования уравнений системы (1)-(4) заключается в том, что при решении системы дифференциальных уравнений (1) относительно концентраций необходимо определять величину И из алгебраического уравнения (3). Существенно затрудняет решение и тот факт, что вместо обычных краевых условий приходится использовать интегральные условия (4). Тем не менее для сравнительно небольших значений J указанная задача без труда решается численно и дает гауссовы распределения концентраций. Однако для больших значений плотности электрического тока

J перед функциями С в уравнениях (1) появляется

к

большой параметр. Это приводит к сильным изменениям производных при слабых изменениях С , что

к

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

Преобразование системы

Оказывается, что задача (1)-(4) допускает замену переменных, позволяющую свести ее к обычной краевой задаче.

Введем в рассмотрение новую функцию ^, такую, что И = кк - е^, где к = 10 - квадратный корень из

w

ионного произведения воды. Кроме того, примем обозначения: ^ = 1Л1/к(к) - К(к) /к2], 5 = 0,5-,/К{к)/К{к) .

к 2 V 1 2 к) к 2

Пусть теперь С"™ - новые неизвестные функции, свя-

к

занн^1е со старыми функциями С соотношениями:

к

С = 2к - . Кроме того, пусть J = 2к - Jnew. Если

к к к к

„new

теперь представить % в виде произведения двух но-

k

вых функций с и ср (у):

new

£ = с ф (у)..

k k k

(5)

где ер (у) = S^ + ch(у-у), то задача (1)-(4) сведется к системе одномерных дифференциальных уравнений: £, (x) из формулы:

чина ( « 0,0257 )) обеспечивал высокую точность решения: малым изменениям функции F (x) соответст-

k

вовали малые изменения величины c .

k

Функции F (x), как следует из уравнений (5), (9),

k

позволяют определить концентрации амфолитов

dc jnew dp c k . J k k

-s—- +-

dx y

dn

— = с ■p ;

dx k k

dx p

= 0 ;

£ = 2k (д + ек(у —у ))• Ь ехр| -^ I.

k w k k k I е k I

Соответственно, функция проводимости смеси у может быть найдена из формулы

N

y = S л

k=1 k

у = — ■ ln

2

( А 2 d р

dx

Р

dp

k

dx

2

ск+/ch ; (6)

N

y = S л k=i k

sh (у-у )

ch (у-у )--k—

k S + ch {у-у )

b expf — F ) + k 4 s k I

N

1 + S С ■ e

у,

k=1

N

1+ S С ■ e

\-1 ^

- у

k

k=1

n (0) = 0 ;

k

n (l ) = m ; k k

k = 1,2,...,N.

(7)

Обратим внимание, что в систему введена новая вспомогательная функция

x

n ( x) = j a ■p (x)dx,

k q k k

представляющая собой количество концентрации k -го амфолита на отрезке [0, x]. Рассмотрение новых функций п , k = , позволило избавиться от РУнге-Кутта (полученные при этом зииенм п (1),

+2^ ек (у — у ), а функция рН - из соотношения:

w k

рН = — ^ • ехр(у)) .

w

Неизвестные функции F (х) определялись из со-

k

ответствующей краевой задачи, для решения которой были разработаны 2 численных алгоритма.

Первый основывался на модифицированных методах Рунге-Кутта и Ньютона и обеспечивал эффективное решение начально-краевой задачи при заданной плотности тока J . Он состоял из 3 этапов. На 1-м этапе производилось задание начальных приближений F (0) = х (k = 1,2,...,N), а также осуществля-k k

лось решение задачи (5)-(7) на отрезке [0,1] методом

k

главной трудности исходной задачи - интегральных условий (4). Укажем также, что при записи системы (5)-(7) использованы очевидные соотношения:

k k а -а = 1 2

Чу-у) 1 dpк{у)

S+ ch {у-W^) Pfo} dx

(8)

Численный эксперимент

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

Во избежание получения отрицательных решений

неизвестная функция е была представлена в виде

k

экспоненциальной зависимости:

с = Ь • ехр| -• F (х) I, k = 1,2,..^ , (9)

k k ^е k )

где Ь - постоянный параметр (в приведенных расче-

k

тах был принят за 1). Параметр 1/е (е - малая вели-

вообще говоря, не удовлетворяют условиям (7)), на 2-м задача F (x,x ) = 0 решалась модифицированным методом Рунге-Кутта, на 3-м выполнялось

решение задачи (5), (6) для точных значений x , по-

k

лученных на предыдущем этапе; одновременно вычислялись £ , pH и y, а также строились графики

k

соответствующих функций.

Второй алгоритм, основанный на так называемом методе движения по параметру, позволял выполнять расчеты в широком диапазоне плотностей тока без существенного накопления вычислительной погрешности. Он также состоял из 3 этапов. На 1-м этапе задавались параметры ячейки и характеристики амфолитов, на основании заданных величин рассчитывались постоянные значения p , S , р , / и асимптотически -

k k o

начальные приближения x . На 2-м к начальной плотности тока J применялся алгоритм 1, строились гран

фики £ , pH и y, на 3-м этапе выполнялось цикличе-

k

ское движение по плотности тока J = J +AJ с при-

н

менением на каждом шаге алгоритма 1.

Программа, реализующая алгоритмы, написана на языке Turbo Pascal 7.0 с использованием стандартного модуля Graph. Результат её работы для каждой

1

\

конкретной системы ИЭФ - серия рисунков, позволяющих исследовать зависимость профилей концентраций амфолитов, рH и проводимости смеси у от плотности тока в широком диапазоне ее изменения.

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

Проведенные расчеты показывают, что при достижении некоторой критической плотности тока система входит в «аномальный» режим: профили концентраций амфолитов теряют сходство с гауссовским распределением. Вначале максимумы на них трансформируются в «плато», а затем сами профили приобретают вид прямоугольников, вплотную примыкающих друг к другу; градиент рH при этом имеет ступенчатый вид. Выявленные тенденции подтверждаются результатами, полученными в [2, 3].

Поскольку результаты получены численно, возникает закономерный вопрос: не является ли «аномальный» режим, наблюдаемый на графиках, результатом накопления вычислительной погрешности (например, в методе Рунге-Кутта) либо иных вычислительных ошибок, приводящих к неверному решению? Для проверки адекватности модели использованы методы асимптотический тестирования [6] и касательных [7]. Результаты, полученные с помощью обоих методов, позволяют сделать вывод об адекватности построенной расчетной модели исходной краевой задаче.

Построение асимптотического решения задачи при высоких плотностях тока

Построение асимптотического решения осуществлено в несколько этапов, потребовавших длительных и скрупулезных расчетов. Решение включало в себя: 1) исключение из системы переменной х и переход к

с как функциям новой переменной ^ ; 2) разложе-

к

ние функций с в ряд по степеням малого параметра

к

к ; 3) получение системы нормальных дифференциальных уравнений для слагаемого нулевого порядка; 4) редукцию системы к одному уравнению в полных дифференциалах с помощью последовательности трехуровневых замен; 5) нахождение решения из условия интегрируемости, возврат к переменной х .

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

п+1 п+1

а —а

1_—1_

a = -a -

n 0

a = a n+1 О

n n

a -a 1 -1

n n \ i n+1 n+1

a -a \-\ a -a 1 -1J i 1 -1

1 N

a = - У m

0 lk=1 k

(10)

где I - длина ЭК; т - общие количества амфолитов,

к

к = 1,2,...Н; разности степеней диссоциации а^ —

определяются уравнениями (8).

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

Т еу(ук)-Я(^к)-ф(^кУ(хк — 2к)

•1 2 Лкг =--

2 (0,5 • Ek-<p(zk ) + 4k (ч )-Ф k

где

E = exp k ^

x - z

kk

2 ^

2a

m

zk = xk + 0,5 • hk +1 a = * 0 ' ' -j2zS

n n \ i n+1 n+1

a -a \-\ a -a

. 1 -1J I 1 -1

Результаты расчетов и их интерпретация

Численная реализация задачи. Рассмотрим несколько примеров расчетов, осуществленных с помощью программы. Расчеты проводились в предположениях: l = 2 дм; r = 0,2 дм; T = 298 К. Плотность тока измерялась в А/дм.2

Пример 1 (рис. 1). Рассмотрена система из 8 абстрактных амфолитов; значения изоэлектрических точек pI заполняют интервал от 4,0 до 7,5 с постоянным шагом ApI = 0,5; разброс значений констант диссоциации (k) (k)

ApK = 2 (рКХ2 = pI ± ApK). Исходные количества

амфолитов одинаковы - M = 0,1 моль. Из графиков

k

следует, что при низких (J = 0,0016 + 0,0064) и средних (J = 0,0128 + 0,0352) плотностях тока профили имеют вид, сходный со стандартным гауссовским распределением. При J = 0,0352 на крайних профилях появляются плато. При высоких плотностях (J = 0,1504 + 0,3552) плато четко видны на профилях всех амфолитов, т.е. система выходит в «аномальный» режим, который не может быть описан гауссовским распределением. В «аномальном» режиме имеет место практически полное расслоение ам-фолитов на прямоугольные области, а также ступенчатый вид градиентов рH и проводимости смеси у .

Пример 2 (рис. 2). В расчетах были использованы характеристики амфолитов-носителей, приведенные в [1]. Произвольным образом выбраны 5 амфолитов с pH < 7 : Asp, М-АБК, а - Asp-His, Tyr-Tyr, IsoGln. Как следует из таблицы, имеет место существенная неравномерность распределения амфолитов по изоэлектриче-ским точкам и константам диссоциации. Исходные количества амфолитов одинаковы - м = 0,1 моль. Рисунок

k

2 показывает, что в процессе эксперимента графики оставались асимметричными, ни один из них не имел вида

a

гауссовского распределения. Однако по мере возраста-

Графическое исследование асимптотического ре-

ния тока четко проявилась прежняя тенденция транс- шения. Асимптотические формулы (10) исследованы формации графиков в «плато», а также и приведения путем сравнения их с результатами, полученными рН и проводимости смеси у к ступенчатому виду. численными методами.

Рис. 1. Численное исследование ИЭФ смеси 8 абстрактных амфолитов

Рис. 2. Численное исследование ИЭФ смеси 5 амфолитов с рН < 7

Характеристики амфолитов-носителей

Амфолит Лк) pKi PK2 pI ApK Коэф. подв. x 104

Asp 1,88 3,65 2,77 1,77 2,97

М-АБК 3,12 4,74 3,93 1,62 3,01

а - Asp-His 3,02 6,82 4,92 3,80 2,11

Tyr-Tyr 3,52 7,68 5,60 4,16 1,56

IsoGln 3,81 7,88 5,85 4,07 2,96

His-His 6,80 7,80 7,30 1,00 1,49

His-Gly 6,27 8,57 7,42 2,30 2,40

His 6,00 9,17 7,59 3,17 2,85

ß-Ala=His 6,83 9,51 8,17 2,68 2,30

Tyr-Arg 7,55 9,80 8,68 2,25 1,58

Пример 3 (рис. 3). Для расчетов произвольным образом выбраны 5 амфолитов с pH>7: His-His, His-Gly, His, pt-Ala-His, Tyr-Arg. Исследование показало, что построенная асимптотика (10) достаточно точно отражает динамику расслоения амфолитов при возрастании плотности тока; при этом приближение тем точнее, чем выше плотность тока. Практически полное совпадение расчетного и асимптотического профиля амфолита наблюдалось в «аномальном» режиме, т.е. как только на соответствующем профиле появлялось плато.

Так, например, из рис. 3 следует, что при плотности тока J = 0,0125 асимптотическое решение (10) точно отражает локализацию амфолитов в пространстве, однако имеются расхождения с расчетными профилями. Полное совпадение наблюдается лишь

j' *,»ns

для профиля Tyr-Arg, на котором уже имеется «плато». При J = 0,0285 профиль p-Ala-His утрачивает сходство с гауссовским распределением, на нем появляется «плато» и одновременно наблюдается слияние расчетной и асимптотической кривых, J = 0,0485 соответствует появление «плато» на профиле His, сопровождающееся слиянием его с асимптотикой. При J = 0,093 кривые профилей и асимптотики сливаются всюду, кроме «вершины» профиля His-Gly, на котором «плато» выражено еще недостаточно четко. Наконец при дальнейшем увеличении плотности тока амфолиты расслаиваются на прямоугольники система полностью переходит в «аномальный» режим и наблюдается полное соответствие асимптотики рассматриваемым профилям.

J-0.024

Выводы

Построенная численная модель позволяет подробно исследовать поведение системы ИЭФ в широком диапазоне плотностей тока. Изложим лишь основные выводы, сделанные в процессе исследования реальных и абстрактных систем.

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

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

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

ния амфолитов по рК , рК , в результате чего

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

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

Поступила в редакцию_

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

Литература

1. Ригетти П. Изоэлектрическое фокусирование. Теория,

методы и применение. М., 1986. 398 с.

2. Mosher R.A., Thormann W. High-resolution computer simu-

lation of the dynamics of isoelectric focusing using carrier ampholytes: The post-separation stabilizing phase revisited // Electrophoresis. 2007. № 23. P. 1803-1814.

3. Thormann W., Mosher R.A. High-resolution computer simu-

lation of the dynamics of isoelectric focusing using carrier ampholytes: Focusing with concurrent electrophoretic mobilization is an isotachophoretic process. Research Article // Electrophoresis. 2006. № 27. P. 968-983.

4. Бабский В.Г., Жуков М.Ю., Юдович В.И. Математиче-

ская теория электрофореза: Применение к методам фракционирования биополимеров. Киев, 1983. 202 с.

5. Жуков М.Ю. Массоперенос электрическим полем. Рос-

тов н/Д, 2005. 216 с.

6. Sakharova L.V., Vladimirov V.A., Zhukov M.Yu. Anomalous

pH-gradient in Ampholyte Solution // arXiv 2009. URL: http://arxiv.org/PS_cache/arxiv/pdf0902/0902.3758v1.pdf (дата обращения: 21.02.2009).

7. Сахарова Л.В. Асимптотическое решение задачи мате-

матического моделирования ИЭФ в естественных градиентах рН // Математика. Экономика. Образование: c6. тр. XVI Междунар. конф., VI междунар. симп. «Ряды Фурье и их приложения», 27 мая - 3 июня 2008 г. Ростов н/Д, 2008. С. 122-130.

24 мая 2011 г.