MATHEMATICAL SIMULATION OF ELECTROPHORESIS OF LIVE CELLS BY MULTIELECTRODE DEVICE

全文:

Электрофорезом называется электромеханическое явление перемещения частиц дисперсной фазы в жидкой или газообразной среде под действием внешнего электрического поля. Электрофорез может быть применен для концентрирования живых клеток в малой части того объема жидкости, в котором клетки располагались до начала процесса. Такое концентрирование используется, например, для выделения малого количества раковых клеток из крови с целью диагностики ранней стадии заболевания. Работа одного из современных микроскопических устройств такого рода описана в статье [1]. В настоящей работе мы представляем нашу математическую модель основных физических процессов, определяющих работу этого устройства, и результаты вычислительных экспериментов, показывающие пути повышения его эффективности. Живая клетка, находящаяся в жидкости, по которой течет электрический ток, может быть приближенно рассмотрена как проводящий шар, окруженный тонкой мембраной, являющейся изолятором. При достаточно больших частотах внешнего электрического поля, создающего этот ток, мембрана становится проницаемой за счет тока смещения, и клетка ведет себя как хороший проводник, если, конечно, проводимость окружающей жидкости мала. Отметим, что для соблюдения последнего условия приходится предварительно заменять плазму крови, имеющую проводимость выше, чем внутренняя часть клетки, на раствор сахара, который имеет низкую проводимость, но обеспечивает осмотическое равновесие, без которого клетки гибнут. Соответствующая теория изложена в монографии [2], где получено распределение электрического поля вокруг такого объекта. Возмущение внешнего электрического поля соответствует возникновению дипольного момента p , величина которого зависит от параметров клетки, жидкости и частоты поля. В рассматриваемом микрочипе [1] используется такой набор параметров, что этот дипольный момент лишь на несколько процентов меньше того, который возник бы на идеально проводящем шаре: p = 4πΚ3εε0E, (1) где R - радиус шара; ε - диэлектрическая проницаемость жидкости; ε0 - диэлектрическая проницаемость вакуума; E - напряженность внешнего электрического поля в жидкости. Если электрическое поле неоднородно, на диполь действует сила F = ( p grad ) E, которая с учетом (1) может быть записана в виде F = 2πR3εε0 grad E2 . При движении шара в вязкой жидкости на него действует сила трения, которая может быть вычислена по закону Стокса: T = -6πRηV, где V - скорость шара; η - динамическая вязкость. В рассматриваемых явлениях инерцией можно пренебречь, и из равенства нулю суммы рассматриваемых двух сил получаем вектор скорости движения шара V = ■ R^0 3η grad E позволяющий записать уравнение для координат центра шара r в виде dr dt R2 εε, 3η 0 grad E2(r) . (2) Таким образом, математическая модель процесса сводится к построению пространственного распределения напряженности электрического поля и интегрирования системы обыкновенных дифференциальных уравнений (2) независимо для каждой клетки. Краевая задача для электрического поля. При используемой в моделируемом устройстве [1] частоте 1 МГц в жидкости с относительной диэлектрической проницаемостью ε = 78 и проводимостью σ = 0.00176 См/м длина электромагнитной волны λ » 30 м и толщина скин-слоя δ » 10 м. Поскольку размеры устройства в десяток тысяч раз меньше, мы пренебрегаем вихревым электрическим полем. Его напряженность будет оценена ниже. Для бизвихревого поля можно ввести электрический потенциал V такой, что E = - grad V . (3) Он удовлетворяет уравнению электропроводности, которое в однородном проводнике, занимающем область Ω , имеет вид уравнения Лапласа - Δ V = 0 . (4) Стенки сосуда, ограничивающего жидкость, состоят из N участков изолятора Гп , на которых равна нулю нормальная компонента плотности тока, что эквивалентно нулевой производной потенциала в направлении, нормальном к границе: dV і dv I = 0, n = 1, N , (5) и из M участков идеального проводника, т. е. электродов γ m , на которых заданы значения потенциала V = 0, m = 1, M . (6) В достаточно общих предположениях относительно формы области, граничных изоляторов и электродов краевая задача (4)-(6) имеет, причем единственное, решение. Доказательство - аналогично приведенному в работе [3] для задачи с двумя электродами. Γ n 14 Математика, механика, информатика а б Рис. 1. Форма электродов. Радиус внешнего электрода 525 мкм (а). Электрическое поле на последнем этапе концентрирования клеток: жирными отрезками ниже горизонтальной оси показаны электроды кроме двух внешних, тонкие кривые - эквипотенциали с шагом δV = 1 В, жирные - линии тока. Между соседними линиями течет ток δC = 1 мкА (б) Моделируемый микрочип представляет собой заполненный жидкостью цилиндр диаметром 4,5 мм и высотой 6 мм с пластиковой боковой стенкой и стеклянными крышкой и дном. На нижнем стекле напылены электроды, форма которых показаны на рис. 1, а. Радиус внешнего электрода 525 мкм. Если сделать электроды в виде целых концентрических колец, подавая напряжения на них через отверстия в стекле, задача была бы осесимметричной и потенциал V(r, z) не зависел бы от третьей цилиндрической координаты φ. Соответствующие напряженность электрического поля и плотность тока имеют нулевые азимутальные компоненты. Поэтому решение остается справедливым для любого сектора φ1 < φ < φ2, ограниченного изолятором. Например, можно установить изолятор вместо всей нижней на рис. 1, а половины устройства, исключив тем самым не кольцевые участки электродов, что несложно сделать на практике. Поскольку такая осесимметричная задача гораздо проще решается, чем трехмерная, но позволяет увидеть важные свойства микрочипа, мы ограничимся именно осесимметричной моделью. Уравнение электропроводности (4) принимает вид д DV D2V -—(r—)-—т Dr Dr Dz = 0, (7) ÖV, Dr I r^0 = 0. (8) конечности подведенной к проводнику электрической энергии. Работа устройства состоит из последовательности 8 этапов, на каждом из которых в течение 20 с напряжение подается на пару соседних электродов, начиная с внешних. Соответственно, задачу (7)-(8) необходимо решить 8 раз со своими наборами значений Vm . Наш многосеточный вариационно-разностный метод численного решения таких задач изложен в монографии [4]. Более простой, но достаточно эффективный алгоритм изложен в работе [5]. В дальнейшем целесообразно использовать сетки, описанные в работе [6]. Для наглядности представления результатов наряду с эквипотенциалями V(r, z) = const используем линии тока, которые в нашем случае совпадают с линиями уровня C (r, z ) = const токовой функции: C(r, z) = j jz (r ', z)2πr ' dr '. условия (6) ставятся на M = 9 электродах, условия (5) -на остальной части границы. На оси симметрии, r = 0 , как на дополнительном участке границы следует выполнить условие Чтобы его получить, проинтегрируем закон сохранения заряда по малому цилиндру, окружающему ось, в предположении гладкости функций: πr 2σ( Ez (0, z + δz ) - Ez (0, z)) + 2πrδzσ Er (r, z) = 0. При r ^ 0 получаем Ez (0, z + δz) - Ez (0, z) = -2δzEr (r, z) / r , т. е. в некоторых точках этого интервала на оси Ez = œ, если Er (0, z ) Φ 0, так как 1/r ^œ. Такое поле соответствует логарифмической особенности плотности джоулевой диссипации, что противоречит На рис. 1, б показано электрическое поле на последнем этапе концентрирования клеток, когда на центральном электроде потенциал V1 = 10 В, на соседнем V2 = -10 В, и Vn = 0 - на остальных. На рис. 2 тонкими линиями показаны линии уровня E2 для первого этапа, когда V1 =... = V7 = 0, V8 = 10 В, V9 = -10 В. Следует отметить, что вблизи всех краев электродов поле имеет особенность E » 1/p^-œ, где ρ - расстояние от линии, разделяющей электрод и изолятор. Эта линия изображается точкой на наших рисунках в плоскости r, z . Разумеется, реальное поле не растет до бесконечности, хотя это и не противоречит конечности джоулевой диссипации. Возрастание ограничено конечностью толщины электродов. Она намного меньше диаметра клеток 2R = 10 мкм, который мы задаем, следуя статье [1]. Поэтому эти особенности нас пока не интересуют. Получившиеся электрические токи C < 10~4 А. Они порождают азимутальное магнитное поле B(r, z) = μ^(r, z)/(2π) <10-3 Тл. Вихревое электрическое поле удовлетворяет уравнению индукции 15 Вестник СибГАУ. № 4(50). 2013 Edl = —j BdS . 5tJ Поскольку магнитное поле велико лишь в небольшой области вблизи пары электродов, имеющей размер порядка 30 мкм, и частота f = 1 МГц, получаем оценку правой части < 10-3 В. Эта величина действительно пренебрежимо мала по сравнению с интересующими нас потенциалами порядка вольта. Этим оправдано представление (3). Движение клеток. Хотя изначально клетки распределены по жидкости случайным образом, чтобы охарактеризовать эффективность их перемещения в малую область, можно расставить их регулярно с постоянной по объему плотностью, как показано светлыми кружками на рис. 2. С увеличением расстояния от оси клетки расположены чаще, поскольку пропорционально r растет объем, приходящийся на каждую нарисованную клетку. Численное интегрирование системы уравнений (2) дает траекторию центра каждой клетки. Используем метод Эйлера с пересчетом. Полученные траектории показаны жирными линиями на рис. 2. Темными кружками показано положение клеток по окончании первого этапа. Изображена только та часть устройства, в которой перемещение клеток заметно. Слияние траекторий не означает, что клетки сталкиваются, они лишь оказываются на одной окружности. На рис. 3 показано положение клеток по завершении восьмого этапа и проведены траектории за все этапы концентрирования, включая начерченные на предыдущем рис. 2. Тонкие линии - это линии уровня E2 для последнего этапа. Они соответствуют приведенным на рис. 1, б эквипотенциалям. Около всех краев электродов есть особенности поля, которые втягивают клетки. Поэтому клетки, оказавшиеся после некоторого этапа вблизи края любого электрода, так и остаются там на следующих этапах, и значит, не приходят в нужную нам центральную область. Этого можно избежать, если сделать для клеток недоступный приграничный слой толще 3 мкм, как это сделано при расчете траекторий, показанных на рис. 3. Рис. 2. Траектории движения клеток: светлые кружки - начальное положение, темные - после первого этапа; тонкие кривые - линии уровня E2 в логарифмическом масштабе для первого этапа, сплошные - при E2 > (104 В/м)2 и штриховые - при меньших напряженностях. Значения E2 на соседних линиях отличаются в VÏÔ раз Рис. 3. Траектории движения клеток в течение всех восьми этапов: светлые кружки - начальное положение, темные - после восьмого этапа; тонкие линии - линии уровня E2 в логарифмическом масштабе для восьмого этапа, сплошные - при E2 > (104 В/м)2 ; значения E2 на соседних линиях отличаются в VÏÔ раз 16 Математика, механика, информатика Все захваченные клетки оказались вблизи внешнего края второго электрода. Видные на рисунке особенности поля занимают примерно по половине этого электрода. Поэтому только дополнительный уход на 15 мкм влево позволил бы клеткам перескочить на левый край электрода. Удалось собрать 41 клетку из 100 изображенных, т. е. всего 0,35 % клеток из 12 000 имевшихся во всем объеме жидкости, который составляет 95 мм3. Объем эффективного захвата клеток мал соответственно глубине проникновения электрического поля, которая определяется масштабом электродов. Это общий недостаток устройств с неподвижной жидкостью. В проточных устройствах типа [7] клетки приближаются к электродам с потоком жидкости, а сильное поле удерживает их в малой области, не мешая жидкости утекать дальше. Для увеличения доли выделенных из жидкости клеток можно уменьшать ее объем. Расчеты показывают, что уменьшение диаметра сосуда до 2 мм и высоты до 0,5 мм распределения полей в областях, показанных на рис. 1-3, практически не изменяются, и поэтому количество собранных на второй электрод клеток остается прежним. Поскольку полное число клеток в жидкости при той же концентрации уменьшается в 60 раз, доля собранных клеток увеличивается до 20 %. Построенные распределения поля также позволяют найти плотность джоулевой диссипации, которая определяет нагрев жидкости. Сила тока, протекающего через устройство, составляет от 80 мкА на первом этапе до 15 мкА на последнем. За весь процесс выделяется около 0,03 Дж тепловой энергии, которая обеспечивает нагрев жидкости на 0,1 К. Если массу жидкости уменьшить в 60 раз, нагревание составит примерно 5 К, что еще допустимо. Нагревание происходит, в основном, вблизи электродов. Соответствующее тепловое расширение жидкости приводит к возникновению конвекции [8], что должно быть учтено при совершенствовании модели. В некоторой мере конвекция уменьшает роль отмеченного выше недостатка конструкции, так как обеспечивает конвективную доставку клеток к электродам. Таким образом, построены пространственные распределения электрического поля и рассчитаны траектории движения клеток в осесимметричном многоэлектродном устройстве, предназначенном для электрофореза живых клеток. Показано, что эффективность его работы ограничена малостью пространственного масштаба поля и наличием особенностей поля, которые возникают вблизи краев тонких электродов.

作者简介

Valery Denisenko

Institute of Computational Modelling RAS SB

Email: denisen@icm.krasn.ru
Doctor of Phisical and Mathematical Sciences, Professor, Leading researcher 50/44, Academgorodok, Krasnoyarsk, 660036, Russian Federation

Sergey Zamay

Krasnoyarsk State Medical University named after prof. V. F. Voino-Yasenetsky

Email: sergey-zamay@yandex.ru
Candidate of Phisical and Mathematical Sciences, main specialist 1, Partizana Zheleznyaka str., Krasnoyarsk, 660022, Russian Federation

Anna Zamay

Krasnoyarsk State Medical University named after prof. V. F. Voino-Yasenetsky

Email: annazamay@yandex.ru
Candidate of Biological Sciences, associate professor 1, Partizana Zheleznyaka str., Krasnoyarsk, 660022, Russian Federation

参考

补充文件