Метод моментов при расчете параметров каналов в микроразмерных системах

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

ВВЕДЕНИЕ                     Данная работа посвящена анализу и моделиро ванию процессов, происходящих в каналах мик-

В последнее время в аналитической технике все чаще используются микро- и нанотехнологии, позволяющие создать современные высокоэкспрессные системы анализа веществ. В частности, интенсивное развитие получили аналитические системы, известные как " ц -TAS" — microTotal Analytical System или "lab-on-a-chip". Основой таких систем является микроразмерная система — микрофлюидный чип, на котором осуществляются все действия с микроколичеством жидкой пробы, включая ввод и дозирование пробы, перемешивание, смешивание с реагентами, реализацию химических реакций, сепарацию полученного продукта, детектирование, сбор фракций. При изготовлении микроразмерных систем (МС) используются прогрессивные методы микро- и нанотехнологий (методы микроэлектронной техники, методы прецизионной литографии, LIGA-технологии и т. д.). Микрофлюидный чип (МФЧ) представляет собой компактное устройство с разветвленной системой микроканалов, микрососудов, реакторов и прочих функциональных элементов. В процессе анализа каналы и сосуды МФЧ заполняются пробой, реагентом, буфером, полимерами. Управление движением микропотоков жидкости осуществляется гидравлическим или электрокинетическим способом. При этом в каналах МФЧ вещество (проба, компонента смеси, реагент) осуществляет сложное (конвективное и диффузионное) движение, что приводит к эффектам, оказывающим влияние как на движения микропотоков, усложняя процессы управления ими, так и на электрофоретическое разделение пробы, искажая результаты анализа.

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

ГЕОМЕТРИЧЕСКИЕ (КОНСТРУКТИВНЫЕ) ПАРАМЕТРЫ КАНАЛА МФЧ И ФИЗИЧЕСКИЕ ХАРАКТЕРИСТИКИ АНАЛИЗИРУЕМЫХ

ОБЪЕКТОВ

В качестве исходных условий полагаем:

• •    что канал МФЧ является прямолинейным и прямоугольного сечения с полушириной h = = 20 мкм и глубиной, равной h (соотношение ширины к глубине, как 2 к 1);

• •    длина канала L является определяемой величиной, но она многократно превосходит h : L / h >  10 2 -10 3 ;

• •    средние значения конвективной скорости U = 10^–2–10^–1 см/с;

• •    коэффициент диффузии вещества D имеет порядок 10^–8–10^–7 см²/с.

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

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

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

Обозначив аксиальную (маршевую) координату x , перпендикулярную стенке координату у и время t , применим пространственное двумерное уравнение Навье—Стокса, описывающее массоперенос вещества в каждом слое без учета реакции, в форме

— + U ( у ) — = D д t         д x

д ² C д ² C

+ дx2   ду 2

Здесь C = C(x, у, t) — распределение концентрации вещества; U(у) — профиль конвективной скорости; D — коэффициент диффузии вещества. Гра-d C (x ,0, t) n ничные условия: симметрия --------= 0 и не ду д C (x, h, t) проницаемость стенки ---------= 0, а также ес- ду тественное условие недоступности вещества ,        дpC(+” у,t)                     „ в форме --------------= 0, p = 0,1,2,... Послед- дxp нее условие говорит о невозможности ухода анализируемого вещества на бесконечное (очень большое) расстояние от зоны анализа. Начальные условия связаны с характером заполнения канала пробой и с объемом пробы.

Используем нормированные координаты и параметры z = у / h ; u * = U / h ; D * = D / h ². Здесь U — максимальная скорость конвективного движения (на оси канала). Полагаем скоростной профиль и ( z ) = и *(1 - | z| m ); z g [0;1]. Параметр m , характеризующий степень клиновидности профиля, априорно неизвестен. Примем, что m > 2 (случай m = 2 — вариант параболического профиля, известного в хроматографии). Моделирование скоростного профиля и оценка параметра m требуют решения уравнения Пуассона—Больцмана для электростатического потенциала, температуропроводности и профиля конвективной скорости.

Решение этой задачи при упрощающем положении о равенстве температуры по всему каналу будет приведено в разделе "Моделирование скоростного профиля электроосмотического потока".

Для решения уравнения (1) воспользуемся методом моментов, предложенным в работе [1]. Решение аналогичной задачи применительно к хроматографии описано во многих работах, например в [2].

Введем моменты порядка п , определяемые как

+^

^ n ( z , t ) = J x"C ( x , z , t )d x .              (2)

-^

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

Момент р имеет смысл количества вещества, отношение д = р ₁/ р ₀ определяет центр тяжести аналитического сигнала, его дисперсия определяется как о ² = р ₂/ р ₀ - п ². Вычисление третьего начального, а затем и центрального моментов, позволяет оценить асимметрию аналитического сигнала (распределения вещества). Центральный момент четвертого порядка дает оценку, аналогичную коэффициенту эксцесса распределения, что используется при проверке гипотезы о гауссовом характере распределения вещества.

Исходное уравнение для моментов в нормированных параметрах примет вид:

^р д t

-

D * •

^{д 2} Р п

d z ²

= пи * (1 - | z| m ) Р п - 1 + п ( п - 1) D * Р п - 2 ; п = 0,1,2,3, ^

Граничные условия для моментов любого порядка являются граничными условиями второго рода: ^{Э р} " ⁽⁰ - ' ) = 0, ^{Э р} " ^{(U 1} = 0. ^, о z             о z

Начальные условия определяются исходным распределением вещества, которое полагаем соответствующим "пробке" относительной длины 2 А с равномерным распределением вещества по сечению, то есть при x G [ -А ; А ] ^ C ( x , z ,0) = 1/(2 А ). Очевидно, что переход к естественным значениям по концентрации произойдет при учете множителя 2 А С ₀. Подобное нормирование удобно с точки зрения формирования начальных условий для моментов:

р 0 ( z ,0) = 1; Р 1 ( z ,0) = 0;

р ₂( z ,0) = А ²/3; р з ( z ,0) = 0; ^

ОПРЕДЕЛЕНИЕ ЦЕНТРА ТЯЖЕСТИ И ДИСПЕРСИИ ПИКА ИНФОРМАТИВНОГО АНАЛИТИЧЕСКОГО СИГНАЛА

Решение полученных уравнений осуществляется по следующему алгоритму: для n = 0 уравнение является однородным и решается методом разделения Фурье; для n = 1, 2,… правая часть уравнения должна быть разложена по собственным функциям T j ( z ) (они же — решения задачи Штурма—Лиувилля, соответствующие уравнению для Ц ).

дцп    * д2ц

При n = 0 имеем ---= D --— и для функ- д t        дz1

ций T j ( z ) получаем уравнения: T " + X T = 0;

T ' (0) = 0, T ' (1) = 0. Поэтому собственные функции (решения задачи Штурма—Лиувилля) будут T j ( z ) = cos( n j z ); j = 0,1,2,...

Исходя из разложения начального условия в ряд Фурье по косинусам, ц ₀ ( z , t ) = 1. Этот результат означает равномерное распределение вещества по всей ширине канала.

Для последующих этапов решения требуется разложить конвективный скоростной профиль по функциям T j ( z ), то есть фактически применить разложение в ряд Фурье по косинусам.

В дальнейшем будем полагать, что

+^

1 - z m = B ⁰ + ^ B^j cos( n j z ).

j = 1

Коэффициенты Фурье определяются параметром m и приведены в табл. 1.

При больших m коэффициенты вычисляются:

B

m m +1

и

Табл. 1. Коэффициенты ряда Фурье разложения профиля конвективной скорости

m	B 0	B j
2	0.667	4( - 1) j + 1 /( n j )2
3	0.750	6( - 1) j + 1 /( n j )2 - - 12 [1 - ( 1> | . n / T
4	0.800	8( - 1) j + 1 /( n j )2 -- 48( - 1) j + 1 /( n j )4
5	0.833	10( - 1) j + 1 /( n j )2 -- 120( - 1) j + 1 /( П )4 + + 240 |1 - ( - 1) j ] /( n j )6

Система обыкновенных дифференциальных уравнений для функций V j , включающая начальные условия, примет вид:

-V = u* • B °; V ₃(0) = 0 и

- j + D* • ( n j )² V j = u* • B j ; V j (0) = 0.

Отсюда следует, что       V 0 = u * • B ⁰ • t ;

* j

V j = n* / Л2 ¹ - exp^{( -} D* • ^{( n} j ⁾² • t ^{) ]} " D • ⁽ n j )

Среднее значение первого момента и, в нашем случае, центра тяжести пика, равно V 0 .

В случае больших диффузионных времен t оценка центра тяжести будет

BX + ² =

m + 2 (j

( 2 • ( - 1) ^j - ( m + 1) • Bj\_m ) .

+^    * n j u • B

Ц 1 ^{( z} , t ) = u • B • t + >              c^os( n jz ).

7 = 1 D • ( n j ) ²

Таким образом, на основании приведенной рекуррентной формулы можно продолжить последовательности коэффициентов разложения B j для любого целого m (на основе данных табл. 1).

Эта формула верна и для нецелых значений m . Но начальные значения коэффициентов B j нужно вычислять численно (оттабулировав для 0 <  m < 2).

Для случая n = 1 ищем решение в форме

+^

Ц 1 ^{( z} , t ) = V ) ⁽ t ) ⁺ ^ V j ⁽ t )^{cos( n j z} ).

j = 1

При этом формула будет достаточно точна, если в первом слагаемом суммы ( V 1 ) функции exp(...) достаточно малы, например, менее 0.01. Соответствующие диффузионные времена должны ln(100) 0.467 удовлетворять требованию t >>  —2—— = —.

В дальнейшем такие времена будем считать большими.

Возведение в квадрат и усреднение по z (по интервалу [ 0;1 ] ) приводит к оценке

) Д 1 ( t \ = ( и

• B ⁰ • t )² + ¹ £ u ^#2 ' B^j

• ² j = 1 D   • ( П )

  .

  
  

Учитываем, что cos2 (njz )^ = 1/2;

) cos( n jz )( = 0.

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

В этих условиях больших диффузионных времен оценка величины дисперсии пика примет вид линейной функции

)” ²< =[ 2 + Pe² £ B L ] ( j = 1 ^{( n j} ) J

" А ² 3_ ₂^4 B^{j 2 4}

---Pe² 2 -----

^ 3 2    ^X_xn )⁴ J

D * • t +

+

.

Здесь число Пекле определяется как Pe = = u / D . Воспользовавшись формулами для сумм рядов вида Sk

+^

= £ 4г j = 1 ^j

[3], можно оценить соот-

ветствующие

+” D j 2

S 2 = X Bv Wj )

+” B j 2

коэффициенты S 1 = X---- j 1 ( n j )²

и

для различных профилей (парамет ров m) (см. табл. 2).

Табл. 2. Характерные параметры, определяющие дисперсию пика для различных профилей (параметров m )

m	S г102	S 2 ^103
2	1.693	1.693
3	1.389	1.357
4	1.108	1.057
5	0.890	0.829

Анализ выражения (3) позволяет сделать следующие выводы:

• 1)    структура (3) качественно напоминает выражение для тейлоровской дисперсии хроматографического пика. То есть имеется слагаемое, соответствующее исходному "размытию" пробы, ( А ² / 3) а также линейно возрастающие со временем — диффузионное (первое) и конвективнодиффузионное (второе слагаемое в первой скобке);

• 2)    это совпадение неслучайно, учитывая сходство процессов конвективно-диффузионного мас-сопереноса в задачах хроматографического и электрофоретического разделений;

• 3)    по мере возрастания клиновидности профиля ( m ) дисперсия пика убывает, что формально отражено в монотонном уменьшении суммы S 1 с ростом m ;

• 4)    по мере роста параметров А и и происходит возрастание дисперсии пика; зависимость от D не столь однозначная. В частности, в случае D * = и * • S 1 /2 первое слагаемое (3) минимизируется. То есть получаем наличие нетривиального оптимального (конвективно-диффузионносогласованного) режима измерений.

В табл. 3 приведены результаты расчета параметров пика (центра тяжести и дисперсии) для широкого диапазона времен анализа. Расчетные параметры m = 5, D = 0.05, и = 7, А = 5.

Видно, что параметры (центр тяжести и дисперсия пика) в различных точках сечения канала ( z ) не очень существенно отличаются друг от друга. Самым быстрым будет движение пика по оси канала, самым медленным — вблизи стенки. Кроме того, "отставание" центра тяжести пика, формируемого у стенки канала от пика, расположенного по оси, практически не накапливается с течением времени. Благодаря монотонному возрастанию дисперсий пиков уже при временах 100 и более обеспечивается практически полное перекрывание пиков со всех точек сечения. Тем самым для расчетов можно считать достаточным использование средней оценки (3).

Оценка по формуле (3) должна быть занижена. Особенно этот эффект существенен при малых временах. Непосредственный численный расчет значений вторых моментов (квадрата центра тяжести) для 101 равномерно отстоящих с шагом 0.01 значений z с последующим усреднением дает оценки: 72.78 ( t = 10), 161.20 ( t = 20), 249.80 ( t = = 30) и 427.20 ( t = 50). Относительная погрешность численного расчета лежит в пределах 1 %.

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

Табл. 3. Параметры пика (центр тяжести — сверху, дисперсия---снизу курсивом)

t	10	20	30	50	100
z
0.00	61.777,	120.139,	178.472,	295.139,	586.806,
	64.607	152.586	240.835	417.337	858.590
0.40	59.936,	118.278,	176.611,	293.278,	584.944,
	81.081	169.160	257.410	433.911	875.165
0.80	55.061,	113.371,	171.705,	288.371,	580.038,
	68.524	156.522	244.771	421.272	862.525
0.95	53.633,	111.938,	170.271,	286.938,	578.604,
	54.394	142.308	230.550	407.051	848.305
Среднее*	58.333,	116.667,	175.000,	291.667,	583.333,
	72.18	160.40	248.62	425.06	866.16

* Под средней оценкой подразумевается расчет по формуле (3).

В свою очередь, коэффициент асимметрии определяется через центральный момент третьего порядка как A = χ ₃ / σ ³ . Учитывая нормировку анализируемого вещества на единицу ( µ 0 = 1), центральный момент третьего порядка вычисляется как χ ₃ = µ ₃ - 3 µ ₂ µ ₁ + 2 ( µ ₁ ) ³ , где моменты µ i вычисляются по (2). Соответствующее дифференциальное уравнение в частных производных было сведено к системе ОДУ и решено численно для небольших времен t . Результаты приведены в табл. 4.

В целом по мере увеличения t , D ^* и ∆ коэффициент асимметрии убывает до значений, порядка 0.01. Рост первых двух параметров свидетельствует об усилении процессов обмена между сечениями канала. Большие длины пробки формируют изначально симметричное распределение вещества, которое будет полностью доминировать на первом этапе движения вещества и компенсировать неравномерность конвективного потока. Аналогичный эффект дос _* тигается уменьшением конвективной скорости u^* . Таким образом, гипотеза о симметричности пика в первом приближении подтверждается.

ЗАСТОЙНЫЕ (ПРИСТЕНОЧНЫЕ) ЗОНЫ ВВОДА ПРОБЫ

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

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

Cα

C ( x , z ,0) = ⁰             exp( - α z ),

2 ∆ 1 - exp( - α )

где α — параметр, характеризующий неравномерность распределения вещества по ширине канала. При α = 0 модель переходит к модели равномерного заполнения проб, и по мере роста α неравномерность исходного заполнения возрастает. При этом при z ≈ 0 доля вещества в соответствующих сечениях превосходит 1, при z ≈ 1 — существенно меньше 1.

Разложение исходной пробки вещества в ряд Фурье позволяет получить решение уравнения (1) для µ 0 как

2 α ²

µ ₀( z , t ) = 1 +              ×

1 - exp( - α )

^+∞ 1 - ( - 1) ^j exp( - α )

× ∑                  exp( - D ⋅ π ² j ² t )cos( π jz ).

j = 1     α ² + ( π j )²

Табл. 4. Коэффициенты асимметрии пика при различных условиях ввода и движения вещества по каналу

∆	* u	D *	t	A
Зависимость от длины пробки
1	10	0.1	3	0.1270
2	10	0.1	3	0.1178
3	”	”	”	0.1048
5	”	”	”	0.0756
10	”	”	”	0.0276
15	”	”	”	0.0105
Зависимость от конвективной скорости
5	3	0.1	3	0.0070
5	5	”	”	0.0261
”	8	”	”	0.0585
”	15	”	”	0.1022
Зависимость от коэффициента диффузии
5	10	0.05	2	0.4482
”	”	0.1	”	0.1129
”	”	0.2	”	–0.0108
”	”	0.3	”	–0.0438
Зависимость от времени анализа
5	10	0.1	1	0.1301
”	”	”	5	0.0259
”	”	”	10	–0.0347

Видно, что по мере роста t количество вещества в каждом сечении z выравнивается, стремясь к единице. Далее решаются уравнения для n = 1,2, и рассчитываются параметры пиков. В табл. 5 приведены данные, характеризующие динамику распределения вещества по сечениям канала для модельных параметров ∆ = 10; u ^* = 10; D ^* = 0.1; t = 1 и t = 3; α = 0, 0.5, 1 и 2.

На основании приведенных данных можно утверждать, что даже при достаточно неравномерном заполнении веществом канала (α = 2) влияние начальной неравномерности заполнения "пробки" мало сказывается при временах t ≥ 3. То есть, главным образом, на положение пика влияют параметры, определяющие размывание пика, а именно D* и u*. К моменту времени t = 3 даже в случае α = 2 происходит существенное выравнивание распределения вещества по сечениям канала до уровня (100±3 %). При этом распределение центров тяжести пиков по сечениям z практически не зависит от начальной структуры пробки вещества; соответствующие дисперсии пиков отличаются в пределах ±10 % для приведенных значений параметра α от 0 до 2. Например, для z = 0 при t = 3 в случае α = 0 центр тяжести и дисперсия пика будут 27.333 и 45.116, а в случае α = 2 — 27.981

Табл. 5. Динамика распределения вещества по сечениям канала

z	α = 0	α = 0.5	α = 1	α = 2
Время t = 1
0	1.0000	1.0194	1.0706	1.2197
0.40	”	1.0049	1.0195	1.0627
0.95	”	0.9828	0.9342	0.7921
Время t = 3
0	1.0000	1.0026	1.0095	1.0299
0.4	”	1.0008	1.0029	1.0092
0.95	”	0.9975	0.9906	0.9705

Табл. 6. Распределение вещества по сечениям канала

z	δ = 0	δ = 0.1	δ = 0.2
Время t = 2
0	1.0000	1.0303	1.0648
0.30	”	1.0179	1.0382
0.80	”	0.9754	0.9474
0.95	”	0.9699	0.9357
Время t = 5
0	1.0000	1.0016	1.0034
0.30	”	1.0009	1.0020
0.80	”	0.9987	0.9973
0.95	”	0.9985	0.9967

и 47.819 соответственно. При этом колебания дисперсии пика для различных сечений z сильнее в случае α = 0. Можно предполагать, что при дальнейшем росте времени анализа t эффект влияния изначального неравномерного распределения вещества ("экспоненциальные застойные зоны") еще сильнее нивелируется.

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

Введем "барьерную застойную зону". В рамках этой модели полагаем начальное распределение вещества в "пробке" как

C ( x , z ,0) =

^C 0    1

2 ∆ 1 - δ ^,

0 ≤ z ≤ 1 - δ ;

C ( x , z ,0) = 0, z > 1 - δ .

Таким образом, устанавливаем "барьер" для закачивания вещества в "пробку", расположенный на относительной глубине канала δ, отсчитываемой от стенки. Этот слой исходно не содержит вещества. Разлагая исходный профиль концентраций в ряд Фурье и решая уравнение при n = 0, по- лучим распределение вещества по сечению канала в форме

Ц 0 ( z , t ) = 1 + -—-х

1 - 8

vi (-1) ^{j + 1} sin( n j 8 )      .          ₂ .₂ .      . . .

^ХУ --- exp( ^- D , ⁿ j t )cos( ⁿ jz ).

П

Как и для предыдущей модели, очевидна тенденция выравнивания вещества по сечениям с ростом t . Результаты расчетов динамики распределения вещества для различных характеристик барьера 8 = 0.1, 0.2 и 0 (отсутствие барьера) и при заданных и = 7, D = 0.1, А = 5, t = 2 и 5 сведены в табл. 6.

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

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

РАЗРЕШЕНИЕ ДВУХ ПИКОВ

Рассмотрим два вещества_* , характеризующихся конвективной скоростью u ^* 1,2 и коэффициентом диффузии D ^* 1,2 . Полагаем, что сигнал каждой анализируемой компоненты смеси является гауссовской кривой, т. к. в первом приближении была подтверждена гипотеза малой асимметрии. То есть форма i -го пика имеет вид

Ji ⁽ x ) =

( exp

V

i = 1,2.    (4)

равном (соизмеримом) содержании обеих компонент смеси. То есть суммарная мощность информативного сигнала от каждой из компонент будет одинакова.

Сформируем суммарный сигнал f sum ( x ) = f i ( x ) + f 2 ( x ) и определим координаты точек экстремумов. Необходимые требования по адекватному оцениванию сигнала и достоверному разделению пиков можно сформулировать как а) x ux = m 1,2 , б) f ( x min ) <<  f ( X tT )- Практически требуется, чтобы разница между уровнями минимального и максимального сигналов превосходила двойную амплитуду шума.

Можно показать, что для четкого разделения двух соседних равномощных пиков (при отношении сигнал/шум, например, 20:1 по амплитуде), требуется выполнение соотношения |^m1 - ^m 2| >  2-42.0 . В этом случае однозначно идентифицируются как оба максимума суммарного сигнала (практически в точках m 1,2 ), так и минимум сигнала между пиками. То есть оба пика практически достоверно распознаются. В дальнейшем будем для простоты руководствоваться критерием 3 ° .

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

В работе [4] приведены данные для конвективной скорости и коэффициента диффузии одноцепочечных фрагментов ДНК длиной 100, 200, 300, 400, 600 и 800 базовых пар, которые воспроизведены в табл. 7. Геометрические характеристики канала: глубина и полуширина канала равны h = = 20 мкм. Объем пробы — 4 А h ³ или 32 А пиколитров (пл).

Рассмотрим задачу оценки минимальной сепа-рационой длины L канала и соответствующего времени разделения соседних фрагментов смеси 100-200-300-400 bp, 400-600-800 bp.

Полагаем объем пробы 320 пл ( А =10). В соответствии с (3) можно вычислить дисперсию пика. Среднее положение центра тяжести пика определяется как и • B • t .

Оценим требуемую длину канала L по следующей методике: будем рассчитывать время прохода центра тяжести каждого пика через точку x = L ; временными границами пика будем полагать время, за которое данный компонент проходит точки L ± 1.5 ° . Соответствующие времена пересчитываются через скорость прохождения центром тяжести пика требуемого расстояния 1.5 ° . В табл. 8 представлены расчетные данные по разделению фрагментов 100, 200, 300 и 400 bp для точки x = L .

Табл. 7. Скорости и коэффициенты диффузии различных фрагментов ДНК (из [4] )

ДНК (bp)	v , мм/с	D ∙105, мм2/с	* u	D *∙102
100	0.1238	1.530	6.19	3.825
200	0.0928	0.789	4.64	1.973
300	0.0744	0.486	3.72	1.215
400	0.0624	0.356	3.12	0.890
600	0.0488	0.239	2.44	0.598
800	0.0411	0.168	2.06	0.420

Табл. 8. Оценка разделения пиков 100-200-300-400 bp в каналах с различной рабочей длиной

L	Фрагмент, bp	Центр пика, t c	σ	∆ t	Время выхода пика t c ± 1.5 ∆ t
800	100	155	37.30	7.2	144–166
”	200	207	44.52	11.5	190–224
”	300	258	50.39	16.3	233–283
”	400	308	53.73	20.7	277–339
1000	100	194	41.80	8.1	182–206
”	200	259	49.89	12.9	239–279
”	300	323	56.53	18.2	296–350
”	400	385	60.30	23.2	350–420

Табл. 9. Характеристики разделения фрагментов 400-600-800 bp в каналах различной длины

L	Фрагмент, bp	Центр пика t c	σ	∆ t	Время пика t c ± 1.5 ∆ t
1100	400	423	63.3	24.3	386–460
”	600	541	68.0	33.4	491–591
”	800	641	74.1	43.0	576–706
1500	400	577	74.2	28.5	534–620
”	600	738	79.8	39.2	679–797
”	800	874	87.2	50.8	798–950

Итак, при рабочих длинах канала более 1000 (2 см) осуществляется хорошо выявляемое разделение всех 4 составляющих анализируемой смеси. При этом полное время анализа составит около 7 мин.

Аналогичную задачу решим применительно к разделению фрагментов 400-600-800 bp. Результаты приведены в табл. 9. Исходя из данных табл. 9, требуемая сепарационная длина канала для разделения фрагментов 400-600-800 bp — 1500 (или 3 см). При этом время анализа составит 16 мин.

В оценочную формулу (3) входят величины, определяемые коэффициентом m , характеризующим клиновидность профиля. Для приведенных выше расчетов данных табл. 8 и 9 взято среднее значение этого коэффициента m = 5. В принципе, скоростной профиль может обладать еще большей кли-новидностью (например, m = 6, 7 и более). В этом случае разрешение увеличивается, во-первых, за счет возрастания средней конвективной скорости (и абсолютной величины разности конвективных скоростей различных фрагментов), во-вторых, за счет уменьшения дисперсии (через коэффициенты S 1 и S 2 , которые убывают с ростом m ).

Приведем сравнительный расчет для случая m = 6. Благодаря большему значению коэффициента Фурье B ₀ = 6/7 (вместо 5/6) и меньшим значениям остальных коэффициентов B j и соответственно сумм S 1 и S 2 (равных 0.726^.10^–2 и 0.662^.10^–3) увеличивается абсолютная величина разницы скоростей движения центров тяжести различных пиков при уменьшении их дисперсии. Для случая разделения фрагментов _х 100-200-300-400 bp для L = 800 имеем временные положения пика (с учетом размаха ± 1.5 σ ): пик 100 (151 ± 10), пик 200 (201 ± 15), пик 300 (251 ± 21) и пик 400 (299 ± 27). Полное разделение наблюдается для пиков 300 и 400 как критических для общего разделения. В данном случае рабочая длина канала должна быть не меньше 1.6 см. При этом время разделения смеси сокращается примерно до 330 с. Аналогичные данные по разделению фрагментов 400600-800 bp имеют вид (для L = 1250): пик 400 (467 ± 34), пик 600 (598 ± 46) и пик 800 (708 ± 60). Отсутствует перекрывание пиков 600 и 800; рабочая длина канала сокращается до 2.5 см, и время анализа не превосходит 13 мин.

Моделированием пограничного (пристеночного) слоя на основе системы уравнений Пуассона— Больцмана позволит оценить параметр клиновид-ности профиля m и тем самым более точно спрогнозировать режим разделения.

МОДЕЛИРОВАНИЕ СКОРОСТНОГО ПРОФИЛЯ ЭЛЕКТРООСМОТИЧЕСКОГО ПОТОКА

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

Помимо введенного ранее нормированного параметра сечения z , отсчитываемого от оси канала полуширины h ( z ∈ [0;1]), проведем нормировку электростатического потенциала на величину дзета-потенциала ζ 0 , то есть ζ ∈ [0;1] .

Решаем систему уравнений Пуасссона— Больцмана при следующих значимых допущениях.

• 1)    Температура в каждом сечении канала одинакова T = 300 K (основываемся на микроразмерах сечения).

• 2)    Полагаем достаточно малой концентрацию буфера, а именно концентрация С 0 находится в диапазоне 10^–6–10^–4 М.

• 3)    Относительная диэлектрическая проницаемость ε = 80 (воде соответствует значение 81) [5].

• 4)    Величина дзета-потенциала ζ 0 не превосходит 100 мВ (лежит в пределах 10–100 мВ). [2]

В этой ситуации коэффициент β = = ( F ζ ₀)/( RT ) <<  1 и sh ( βζ ) ≈ βζ ; F и R обозначают соответственно постоянную Фарадея и универсальную газовую постоянную.

Уравнение Пуассона—Больцмана, описывающее распределение потенциала по сечению канала, линеаризуется (5) и решается аналитически

• ^{d2 ζ} ₂ = λ ² ζ . (5) d z

Граничные условия: симметрии dζ dz

z = 0

= 0

и электростатического баланса ζ I z =1 = 1 , так как на стенке потенциал совпадает с дзета-потенциалом.

Параметр λ определяется, исходя из выражения

λ ² =

F ² C ₀ h ² ε ₀ ε RT ^,

где ε 0 — электрическая постоянная [5];

ε ₀ =     ₂ ≈ 8.842 ⋅ 10 ^{- 12} ф I м , с — скорость света

4 π с ²

в вакууме.

Распределение относительной величины потенциала описывается выражением ζ = ch(λz) / ch(λ) . Поскольку уравнение для скорости электроосмотического потока (скоростного d2u 2EFh2C профиля) имеет вид     =        0 sh(βζ) с гра-

• d z ²       η

ничными условиями симметрии профиля на оси канала и прилипания к стенке канала ( u = 0), то это уравнение линеаризуется и аналитически решается при тех же допущениях, что и уравнение Пуассона—Больцмана. Здесь Е — напряженность электрического поля, η — динамический коэффициент вязкости.

Основные закономерности поведения решения:

• а)    по мере роста параметра λ степень неравномерности распределения потенциала по сечению канала возрастает;

• б)    поскольку вторая производная (радиус кривизны) скоростного профиля пропорциональна потенциалу, неравномерность распределения потенциала приведет к усилению неравномерности скорости по сечению канала;

• в)    аппроксимировав профиль выражением u max ^{(1 -} I z I ^m ) , где m — показатель, характеризующий клиновидность профиля, можно связать m c заданными условиями, а именно дзета-потенциалом ζ 0 , шириной канала 2 h , напряженностью электрического поля Е , концентрацией буфера С 0 , коэффициентом динамической вязкости η . Для этого нужно решить задачу оптимизации по критерию минимума среднего квадратичного отклонения профилей друг от друга.

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

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

Заданы исходные данные:

η = 10^–5 см²/с, ζ 0 = 10 мВ, Е = 10⁴ В/м, h = = 20 мкм. Численный расчет оптимального по критерию минимума расхождения профилей коэффициента m представлен в табл. 10.

Влияние ширины канала иллюстрируется данными табл. 11 при значениях величин: η = = 10^–5 см²/с, ζ 0 = = 10 мВ, Е = 10⁴ В/м, С 0 = = 1.5∙10^–5 М.

Табл. 10. Оценка параметра клиновидно-сти профиля m для буферных растворов различной концентрации

С 0 ∙105	m
0.1	2.2
0.3	2.7
1	4.2
2	5.8
3	7.0
5	>7

Табл. 11. Оценка параметра клиновидности профиля m для каналов различной ширины

h , мкм	m
5	2.2
10	2.9
15	3.9
20	5.1
25	6.3
30	>7

Таким образом, сделанные ранее при анализе разделения фрагментов смесей допущения относительно возможных значений параметра m обоснованы. Можно подобрать условия, при которых степень клиновидности профиля будет характеризоваться значением m от 2 до 7 и более.

Зависимость параметра клиновидности m от концентрации буфера С ₀ и от полуширины канала h

цесс массопереноса вещества, что полностью подтверждает априорные качественные предположения, в частности уменьшение размывания вещества по мере роста клиновидности (параметра m ).

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

• 4.    Определен численными методами конвективный скоростной профиль и соответствующий ему коэффициент m на основе решения задачи оптимизации по критерию минимального отклонения полученного скоростного профиля от модельного u ( z ) = u ⋅ (1 - z ^m ); z ∈ [0;1] .

Работа выполнена при финансовой поддержке, осуществляемой в рамках научной программы Санкт-Петербургского научного центра РАН "Аналитические приборы на основе микрофлюид-ных технологий" (раздел 2, проект 2, 2003 г.).

В рассмотренной модели предполагалось прямоугольное сечение канала, аппроксимируемое совокупностью слоев одной ширины 2 h . Переход к пространственно двумерной задаче даже для этого случая потребует наложения аналогичного условия электростатического баланса ζ = 1 для дна канала.

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

ВЫВОДЫ

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

• 2.    Исследовано влияние формы симметричного конвективного профиля произвольной клино-видности, характеризуемой параметром m , на про-