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

В. Е. Косенко и др. // Вычислительные технологии. 2009. Т. 14. Вып. 6. С. 19-28.

2. Васильев Е. Н., Никифорова Е. С. Математическое моделирование теплового режима гипертепло-проводного радиатора мощного радиоэлемента // Вестник СибГАУ. 2005. Вып. 6. С. 23-26.

3. ОСТЕРМ СПБ. Термоэлектрические охлаждающие модули [Электронный ресурс]. иКЬ: http://osterm.ru.

References

1. Vasilyev E. N., Derevyanko V. A., Kosenko V. E.

et al. Computational Technologies, 2009. vol. 14, no. 6, p. 19-28.

2. Vasilyev E. N., Nikiforova E. S. Vestnik SibGAU. 2005, no. 6, p. 23-26.

3. OSTERM SPB. Thermoelectric cooling modules. Available at: http://osterm.ru.

© Васильев Е. Н., Деревянко В. А., 2013

УДК 57.08+519.6

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

В. В. Денисенко1, С. С. Замай2, А. С. Замай2

'Институт вычислительного моделирования СО РАН Российская Федерация, 660036, Красноярск, Академгородок 50/44. Е-mail: denisen@icm.krasn.ru 2 Красноярский государственный медицинский университет имени проф. В. Ф. Войно-Ясенецкого Российская Федерация, 660022, Красноярск, ул. Партизана Железняка, 1 E-mail: sergey-zamay@yandex.ru

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

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

MATHEMATICAL SIMULATION OF ELECTROPHORESIS OF LIVE CELLS BY MULTIELECTRODE DEVICE

V. V. Denisenko1, S. S. Zamay2, A. S. Zamay2

'Institute of Computational Modelling RAS SB 50/44, Academgorodok, Krasnoyarsk, 660036, Russian Federation. E-mail: denisen@icm.krasn.ru 2Krasnoyarsk State Medical University named after prof. V. F. Voino-Yasenetsky

1, Partizana Zheleznyaka str., Krasnoyarsk, 660022, Russian Federation E-mail: sergey-zamay@yandex.ru

A mathematical model of the microchip that concentrates live cells in a small volume of liquid by electrophoresis is created. Space distributions of electric field are calculated for a cylindrical conductor with axially symmetrical electrodes. Equipotentials, field lines and lines of constant electric field strength squared are plotted. Cells are simulated as good conducted balls bounded with a low conducting membrane. The membrane has no influence on electric current if the frequency of the electric field is high enough. Then a cell behaves as a good conducting ball. It is polarized and a force appears which pulls the cell into the regions with the strong electric field. Trajectories of the cells are defined for

* Работа выполнена в рамках государственного контракта № 14.512.11.0086 Минобрнауки РФ.

the cell movement that is the result of this force and viscosity of the liquid. Singularities of the electric field distributions are found near the lines which separate boundary isolators and electrodes. The factors which limit effectiveness of the microchip are analyzed.

Keywords: mathematical model, electrophoresis, cell.

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

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

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

Соответствующая теория изложена в монографии [2], где получено распределение электрического поля вокруг такого объекта. Возмущение внешнего электрического поля соответствует возникновению дипольного момента p , величина которого зависит от параметров клетки, жидкости и частоты поля. В рассматриваемом микрочипе [1] используется такой набор параметров, что этот дипольный момент лишь на несколько процентов меньше того, который возник бы на идеально проводящем шаре:

p = 4nR3ee0E, (1)

где R - радиус шара; е - диэлектрическая проницаемость жидкости; е0 - диэлектрическая проницаемость

вакуума; E - напряженность внешнего электрического поля в жидкости.

Если электрическое поле неоднородно, на диполь действует сила

F = (p grad) E, которая с учетом (1) может быть записана в виде

F = 2nR3ss0 grad E2 .

При движении шара в вязкой жидкости на него действует сила трения, которая может быть вычислена по закону Стокса:

Т = -6жЯц V,

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

V = ■

Я1єє0

3П

grad E

позволяющий записать уравнение для координат центра шара г в виде

йГ йі

R2 єє,

3n

0 grad E 2(F) .

(2)

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

Краевая задача для электрического поля. При используемой в моделируемом устройстве [1] частоте 1 МГц в жидкости с относительной диэлектрической проницаемостью е = 78 и проводимостью ст = 0.00176 См/м длина электромагнитной волны X и 30 м и толщина скин-слоя 5 и 10 м. Поскольку размеры устройства в десяток тысяч раз меньше, мы пренебрегаем вихревым электрическим полем. Его напряженность будет оценена ниже. Для бизвихревого поля можно ввести электрический потенциал V такой, что

E = - grad V . (3)

Он удовлетворяет уравнению электропроводности, которое в однородном проводнике, занимающем область Q , имеет вид уравнения Лапласа

- Д V = 0 . (4)

Стенки сосуда, ограничивающего жидкость, состоят из N участков изолятора Гп , на которых равна нулю нормальная компонента плотности тока, что эквивалентно нулевой производной потенциала в направлении, нормальном к границе:

дУі

dv I

= 0, n = 1,

(5)

и из М участков идеального проводника, т. е. электродов у т , на которых заданы значения потенциала

V

4m

= 0, m = 1,

М .

(б)

В достаточно общих предположениях относительно формы области, граничных изоляторов и электродов краевая задача (4)-(6) имеет, причем единственное, решение. Доказательство - аналогично приведенному в работе [3] для задачи с двумя электродами.

г

n

а б

Рис. 1. Форма электродов. Радиус внешнего электрода 525 мкм (а). Электрическое поле на последнем этапе концентрирования клеток: жирными отрезками ниже горизонтальной оси показаны электроды кроме двух внешних, тонкие кривые - эквипотенциали с шагом 8У = 1 В, жирные - линии тока. Между соседними линиями течет ток 8Є = 1 мкА (б)

Моделируемый микрочип представляет собой заполненный жидкостью цилиндр диаметром 4,5 мм и высотой 6 мм с пластиковой боковой стенкой и стеклянными крышкой и дном. На нижнем стекле напылены электроды, форма которых показаны на рис. 1, а. Радиус внешнего электрода 525 мкм.

Если сделать электроды в виде целых концентрических колец, подавая напряжения на них через отверстия в стекле, задача была бы осесимметричной и потенциал V(г, х) не зависел бы от третьей цилиндрической координаты ф. Соответствующие напряженность электрического поля и плотность тока имеют нулевые азимутальные компоненты. Поэтому решение остается справедливым для любого сектора ф1 < ф < ф2, ограниченного изолятором. Например, можно установить изолятор вместо всей нижней на рис. 1, а половины устройства, исключив тем самым не кольцевые участки электродов, что несложно сделать на практике. Поскольку такая осесимметричная задача гораздо проще решается, чем трехмерная, но позволяет увидеть важные свойства микрочипа, мы ограничимся именно осесимметричной моделью. Уравнение электропроводности (4) принимает вид

д дГ д 2V

-—(r—)-—Т

дr дr дг2

= 0,

(7)

д¥,

дr

| r ^0

= 0.

(8)

конечности подведенной к проводнику электрической энергии.

Работа устройства состоит из последовательности 8 этапов, на каждом из которых в течение 20 с напряжение подается на пару соседних электродов, начиная с внешних. Соответственно, задачу (7)-(8) необходимо решить 8 раз со своими наборами значений Vm . Наш многосеточный вариационно-разностный метод численного решения таких задач изложен в монографии [4]. Более простой, но достаточно эффективный алгоритм изложен в работе [5]. В дальнейшем целесообразно использовать сетки, описанные в работе [6].

Для наглядности представления результатов наряду с эквипотенциалями V(r, z) = const используем линии тока, которые в нашем случае совпадают с линиями уровня C (r, z) = const токовой функции:

C(r, z) = I jz (r', z)2nr' dr'.

условия (6) ставятся на М = 9 электродах, условия (5) -на остальной части границы. На оси симметрии, г = 0 , как на дополнительном участке границы следует выполнить условие

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

пг 2ст( Ех (0, х + 8х) - Ех (0, х)) + 2лг5хст Ег (г, х) = 0.

При г ^ 0 получаем

Ех (0, х + 5х) - Ех (0, х) = -25хЕг (г, х)/ г ,

т. е. в некоторых точках этого интервала на оси Ех = го, если Ег(0,х) Ф 0, так как 1/г ^го. Такое

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

На рис. 1, б показано электрическое поле на последнем этапе концентрирования клеток, когда на центральном электроде потенциал V1 = 10 В, на соседнем V2 = -10 В, и Vn = 0 - на остальных.

На рис. 2 тонкими линиями показаны линии уровня Е2 для первого этапа, когда ¥1 =... = V-! = 0, У8 = 10 В, V9 = -10 В. Следует отметить, что вблизи всех краев электродов поле имеет особенность Е и 1/ р^го, где р - расстояние от линии, разделяющей электрод и изолятор. Эта линия изображается точкой на наших рисунках в плоскости г, х . Разумеется, реальное поле не растет до бесконечности, хотя это и не противоречит конечности джоулевой диссипации. Возрастание ограничено конечностью толщины электродов. Она намного меньше диаметра клеток 2Я = 10 мкм, который мы задаем, следуя статье [1]. Поэтому эти особенности нас пока не интересуют.

Получившиеся электрические токи С < 10~4 А. Они порождают азимутальное магнитное поле В(г, х) = ц0С(г, х)/(2пг) <10-3 Тл. Вихревое электрическое поле удовлетворяет уравнению индукции

т = —Г Всё. д/-1

Поскольку магнитное поле велико лишь в небольшой области вблизи пары электродов, имеющей размер порядка 30 мкм, и частота / = 1 МГц, получаем

оценку правой части < 10-3 В. Эта величина действительно пренебрежимо мала по сравнению с интересующими нас потенциалами порядка вольта. Этим оправдано представление (3).

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

Численное интегрирование системы уравнений (2) дает траекторию центра каждой клетки. Используем метод Эйлера с пересчетом. Полученные траектории

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

На рис. 3 показано положение клеток по завершении восьмого этапа и проведены траектории за все этапы концентрирования, включая начерченные на предыдущем рис. 2. Тонкие линии - это линии

уровня Е2 для последнего этапа. Они соответствуют приведенным на рис. 1, б эквипотенциалям.

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

Рис. 2. Траектории движения клеток: светлые кружки - начальное положение, темные - после первого этапа; тонкие кривые - линии уровня Е2 в логарифмическом масштабе для первого этапа, сплошные - при Е2 > (104 В/м)2 и штриховые - при меньших напряженностях. Значения Е2 на соседних линиях отличаются в л/Тс раз

Рис. 3. Траектории движения клеток в течение всех восьми этапов: светлые кружки - начальное положение, темные - после восьмого этапа; тонкие линии - линии уровня Е2 в логарифмическом масштабе для восьмого этапа, сплошные - при Е2 > (104 В/м)2; значения Е2 на соседних линиях отличаются в -\/Ї0 раз

Все захваченные клетки оказались вблизи внешнего края второго электрода. Видные на рисунке особенности поля занимают примерно по половине этого электрода. Поэтому только дополнительный уход на 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], что должно быть учтено при совершенствовании модели. В некоторой мере конвекция уменьшает роль отмеченного выше недостатка конструкции, так как обеспечивает конвективную доставку клеток к электродам.

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

Библиографические ссылки

1. Huang C.-T., Amstislavskaya T. G., Chen G.-H. at al. Selectively concentrating cervical carcinoma cells from red blood cells utilizing dielectrophoresis with circular ITO electrodes in stepping electric fields // J. of medical and biological engineering. 2012. Vol. 33(1). P. 51-58.

2. Jones T. B. Electromechanics of particles. New York: Cambridge University Press, 1995.

3. Денисенко В. В. Энергетический метод решения трехмерной смешанной краевой задачи для эллипти-

ческого уравнения с несимметричной матрицей коэффициентов. 1981. 27 с. Деп. ВИНИТИ. № 5647-81.

4. Денисенко В. В. Энергетические методы для эллиптических уравнений с несимметричными коэффициентами. Новосибирск : Изд-во СО РАН, 1995. 204 с.

5. Shaidurov V. V., Tobiska L. The convergence of the cascadic conjugate-gradient method applied to elliptic problems in domains with re-entrant corners // Mathematics of Computation. 1999. Vol. 69, № 230. P. 501-520.

6. Киреев И. В., Ищенко А. В. Фрактальный алгоритм построения двумерных вложенных сеток. // Вестник СибГАУ. 2011. Вып. 2 (35). С. 20-24.

7. Gao J., Yin X.-F., Fang Z.-L. Integration of single cell injection, cell lysis, separation and detection of intracellular constituents on a microfluidic chip // Miniaturisation for chemistry, biology & bioengineering. 2004. Vol. 4. P. 47-52.

8. Andreev V. K., Bekezhanova V. B. Instability of the equilibrium state of liquid with a free surface in the presence of volume heat sources // Fluid Dynamics. 2008. Vol. 43, no. 2. P. 176-184.

References

1. Huang C.-T., Amstislavskaya T. G., Chen G.-H., Chang H.-H., Chen Y.-H., Jen C.-P. Selectively concentrating cervical carcinoma cells from red blood cells utilizing dielectrophoresis with circular ITO electrodes in stepping electric fields. Journal of medical and biological engineering. 2012, vol. 33(1), p. 51-58.

2. Jones T. B. Electromechanics of particles. New York, Cambridge University Press, 1995.

3. Denisenko V. V. Energeticheskij metod resheniya trekhmernoi smeshannoi kraevoi zadachi dlya ellip-ticheskogo uravneniya s nesimmetrichnoi matritsei koeffi-cientov (Energy method for three dimensional boundary value problem of mixed type for an elliptical equation with nonsymmetrical coefficient matrix). Dep. VINITI. 1981, № 5647-81, 27 p.

4. Denisenko V. V. Energetichesky metod dlya ellip-ticheskikh uravnenii s nesimmetrichnymi koefficientami (Energy methods for elliptic equations with nonsymmetri-cal coefficients). Novosibirsk, Publ. house of the Russian Academy of Sciences Siberian Branch, 1995, 204 p.

5. Shaidurov V. V., Tobiska L. The convergence of the cascadic conjugate-gradient method applied to elliptic problems in domains with re-entrant corners. Mathematics of Computation. 1999, vol. 69, № 230, p. 501-520.

6. Kireev I. V., Ishchenko A. V. Vestnik SibGAU. 2011, no. 2(35), р. 20-24.

7. Gao J., Yin X.-F., and Fang Z.-L. Integration of single cell injection, cell lysis, separation and detection of intracellular constituents on a microfluidic chip. Miniaturisation for chemistry, biology & bioengineering. 2004, vol. 4, p. 47-52.

8. Andreev V. K., Bekezhanova V. B. Instability of the equilibrium state of liquid with a free surface in the presence of volume heat sources. Fluid Dynamics, 2008, vol. 43, no. 2, p. 176-184.

© Денисенко В. В., Замай С. С., Замай А. С., 2013