CC BY
50
Том 150, кн. 2
Естественные науки
2008
УДК 573.4+001.8
СРАВНЕНИЕ МЕТОДОВ РАСЧЕТА ЭЛЕКТРОСТАТИЧЕСКИХ ВЗАИМОДЕЙСТВИЙ В АЛГОРИТМАХ ДЛЯ ГРАФИЧЕСКИХ АКСЕЛЕРАТОРОВ
Н.И. Акберова, Д.И. Алишева, Е.Д. Изотова, Д.С. Тарасов
Аннотация
Целью данной работы было сравнение различных схем расчетов дальних электростатических взаимодействий на примере средней биологической системы для вычислений на графических акселераторах. Парные методы сравнивались с методами, используемыми для расчета на CPU: метод PME и метод с использованием радиуса обрезания потенциала. За эталонный метод был принят PME. В качестве модельной системы был взят мутантный белок 4SER. Было проведено моделирование молекулярной динамики с различными схемами расчета электростатических взаимодействий. Оценивались структурные параметры (среднее квадратичное отклонение, радиус гирации, поверхность доступная растворителю) для целой системы, площадь полости активного центра. Наилучшие результаты для мутанта 4SER показали парный метод с использованием радиуса обрезания потенциала 8 Â и переключающей функцией. Для активного центра наилучшие результаты показал парный метод расчета электростатических взаимодействий, основанный на суммации Вульфа с радиусом обрезания потенциала 9 Â.
Ключевые слова: графические акселераторы, парные методы расчета электростатических взаимодействий, катепсин L.
Введение
В последнее время интенсивно изучается проблема учета электростатических взаимодействий в различных методах молекулярной динамики для моделирования биологических систем [1-3]. Известно, что электростатические взаимодействия вносят существенный вклад даже на больших расстояниях в динамические и структурные свойства моделируемой системы. Однако время расчета электростатические взаимодействия возрастает согласно квадрату с линейным увеличением количества атомов.
Одним из способов для уменьшения времени расчета вычисления дальних электростатических взаимодействий является отрезание данных взаимодействий для частиц, находящихся на расстоянии, большем заданной величины (это расстояние называют радиусом обрезания потенциала).
Однако использование обрезания потенциала для дальних электростатических взаимодействий приводит к нарушению непрерывности потенциала и нарушению сохранения энергии. Чтобы избежать этого, используется специальная переключающая функция [4], которая позволяет плавно привести дальние электростатические взаимодействия к нулю. Но даже при такой схеме расчета электростатических взаимодействий не удается полностью избавиться от раз-
личных проблем, таких, как низкая стабильность моделирования [5], искусственная стабилизация молекул воды на радиусе обрезания, приводящая к высокой вязкости [6], большие колебания заряженных аминокислотных остатков от исходного положения [7].
Было предложено несколько путей для того, чтобы избежать недостатков методов, использующих радиусы обрезания потенциала. Наиболее часто используемые методы - это суммация Эвальда [8, 9], ее более быстрая и эффективная реализация - PME (particle mesh Ewald).
В первом методе частица взаимодействует со всеми остальными частицами в ячейке периодичности и со всеми их образами в бесконечном массиве периодических ячеек. Учитываются вклады в потенциальную энергию пары заряд-заряд, экранирующего и компенсирующего полей, поправки на взаимодействие частицы самой с собой, поправки на окружающую среду (для вакуума). Предполагается, что каждый заряд окружен нейтрализующим распределением зарядов противоположного знака, так называемым экранирующим полем. Сумма точечных зарядов таким образом преобразуется в сумму взаимодействий между зарядами и нейтрализующего распределения. Данное предположение позволяет ускорить сходимость суммы. Однако это искусственно введенный член, поэтому в систему вводится второе распределение заряда - компенсирующее поле, которое отменяет действие первого нейтрализующего распределения.
Во втором методе (PME) полная энергия электростатических взаимодействий получается благодаря вкладам ближних и дальних электростатических взаимодействий частиц, причем дальние электростатические взаимодействия рассчитываются с помощью быстрого преобразования Фурье. Данный метод масштабируется как NlgN.
В настоящее время моделирование систем, представляющих практический интерес (размерами сотни тысяч атомов), требует использования дорогостоящих высокопроизводительных вычислительных систем. В связи с этим рядом авторов [10, 11] предлагалось использование графических акселераторов (GPU) для расчета парных взаимодействий. Для графических акселераторов разработаны эффективные алгоритмы парных взаимодействий, но в настоящее время отсутствуют эффективные алгоритмы для расчетов методом PME. В связи с этим представляет интерес предложенный недавно метод расчета электростатических взаимодействий путем попарной суммации [12]. Исходя из методов суммации [13, 14], для расчета дальних электростатических взаимодействий авторами предложена следующая формула:
vdsf (r ) = qq
erfc(ar ) erfc(aRc )
r R
(erfc(aRc) + 2a exp(-a2Rc2)л
Я2 n1/2 R„
(r - Rc )
Rc (1)
Схемы с использованием данной формулы с определенными параметрами дают практически идентичные результаты методу суммации Эвальда для исследованных систем, таких, как чистая вода, раствор NaCl в воде с низкой и высокой ионной силой, сфера аргона радиусом 6 Â в воде.
Однако до сих пор таких исследований для больших белков не проводилось, хотя по понятным причинам они представляют большой интерес. Было опубликовано несколько работ, в которых проводилось сравнение метода PME
и методов, использующих различные схемы радиусов обрезания потенциала. Однако эти исследования проводились на примере чистой воды [6, 15, 16], ионов [17], пептидов [18-20], небольших белков [21, 22] и для расчетов на CPU (центральное процессорное устройство). Поэтому будет интересно сравнить результаты моделирования молекулярной динамики для белка катепсина, полученные с помощью программы GPAMM [23] для графических карточек с различными схемами расчета, с результатами, получаемыми методом PME на CPU для того же белка. Для выявления наилучшего метода расчета и проверки адекватности расчета молекулярной динамики в качестве объекта исследования был выбран 4SER-мутантный белок катепсин L, обладающий кремний-конден-сирующей активностью [24]. Поскольку данный фермент достаточно большого размера (состоит из 220 аминокислот), его структура определена с использованием рентгеноструктурного анализа, при этом известно положение активного центра. В дальнейшем будет изучен механизм катализа ферментативной реакции, осуществляемой данным белком, поэтому он является оптимальным объектом исследования.
Цель данной работы - оценка, какая из схем расчетов дает наилучшие результаты с помощью GPU в сравнении с методами, используемыми для CPU. За эталонный метод был принят метод PME. Результаты этого исследования будут использованы для того, чтобы установить оптимальную схему моделирования молекулярной динамики белка катепсина для изучения поведения его активного центра.
1. Материалы и методы
1.1. Моделируемая система. Кристаллическая структура мутанта силика-теина A катепсина L была получена из Protein Data Bank (PDB индекс 2VHS). Разрешение 1.5 Â. Белок состоит из 4 доменов. Для исследования был взят 1 домен.
Для расчета на CPU все водороды были добавлены для 1 домена в PDB файл программой Gromacs-3.3.1 [25] с образованием системы, состоящей из 3209 атомов. Полученная система была растворена и центрирована в коробке воды размером 7.786 х 6.701 х 6.386 нм в программе Gromacs-3.3.1 с потенциалом amber99 [26]. В результате система состояла из 33359 атомов. Была произведена минимизация энергии в течение 1 пс также в программе Gromacs-3.3.1 с использованием потенциала amber99 и расчет траекторий молекулярной динамики в течение 1 нс с использованием 1) радиуса обрезания потенциала, равного 8 Â и 2) метода PME для дальних электростатических взаимодействий.
Для расчета на GPU все водороды были добавлены для 1 домена в PDB файл с помощью программы HyperChem 7.0 [27] с образованием системы, содержащей 3206 атомов. В программе VMD 1.8.6 [28] полученная система была растворена в коробке воды размером 80 х 70 х 65 Â - система содержала 35015 атомов. В программе GPAMM [23] на GPU была произведена минимизация энергии системы в течение 1 пс и расчет траекторий молекулярной динамики в течение 1 нс при температуре 300 К с использованием 1) радиуса обрезания потенциала, равного 8 Â, и переключающей функцией, 2) аналогично, но с из-
мененными, как в amber [26], 1-4 взаимодействиями, 3) алгоритма, модифицированного по формуле (1), с радиусом обрезания потенциала 8 Â, 4) аналогично, но с радиусом обрезания потенциала 9 Â.
1.2. Обработка результатов. В работе для оценки стабильности моделирования рассчитывались следующие параметры: среднее квадратичное отклонение (RMSD), радиус гирации, площадь полярной поверхности, площадь поверхности, доступной растворителю, площадь полости активного центра, расстояния между аминокислотными остатками в активном центре. Среднее квадратичное расстояние рассчитывалось для всего белка и отдельно для активного центра.
Среднее квадратичное отклонение вычислялось с помощью Trajectory RMS Tool в программе VMD 1.8.6. Для расчета радиуса гирации, площади полярной поверхности и площади поверхности, доступной растворителю, использовалась программа VEGA ZZ 2.2.0 [29]. Чтобы рассчитать площадь полости активного центра и расстояний между аминокислотными остатками в активном центре, была написана программа change_field_3amino.py.
2. Результаты и обсуждения
Исследуемый белок относится к классу протеаз. Белок 4SER (рис. 1) в отличие от немутантного имеет серию точечных мутаций и одну замену в нук-леотидной последовательности (173ESTESDNN180 на 173ISNNQ177). По сравнению с исходным ферментом он обладает более высокой каталитической активностью, это связано с увеличением размера кармана активного центра. Было проанализировано несколько параметров, чтобы оценить структурную стабильность в течение моделирования молекулярной динамики для данного мутанта и его активного центра. Рассматривалось также динамическое поведение молекулы. Кроме того, были описаны и проанализированы изменения площади полости активного центра и расстояний между аминокислотными остатками в ходе молекулярной динамики.
Был произведен расчет среднего квадратичного отклонения по результатам каждого из методов вычислений, чтобы оценить стабильность моделирования.
На рис. 2-6 представлены результаты, полученные: А - в программе GPAMM с радиусом обрезания потенциала 8 Â; Б - в программе Gromacs-3.3.1 с использованием радиуса обрезания потенциала 8 Â; В - в программе Gromacs-3.3.1 с использованием метода PME; Г - в программе GPAMM методом с измененными, как в amber, 1-4 взаимодействиями; Д - в программе GPAMM с использованием методики, описанной в [12]; Е - как в Д, но с радиусом обрезания 9 Â.
Как видно на рис. 2, наилучшую стабильность моделирования относительно метода PME при расчете на графических акселераторах показывают методы А и Г. Метод расчета Б на обычных акселераторах с использованием радиуса обрезания 8 Â показывает худшие результаты относительно PME и методов А и Г.
Все методы расчета электростатических взаимодействий показывают расширение структуры относительно начальной структуры (рис. 3). После начального расширения методы Б, В и Д относительно стабильны, и их результаты колеблются в пределах 35.4-35.8, 34.6-34.9, 35.5-35.7 Â соответственно. Интересным
Рис. 1. 48БЯ-мутантный белок. Атомами показаны аминокислотные остатки, составляющие активный центр
3,5
й л
а* 1
5 |
И п
* О
а) §
2,5
1,5
0,5 4
1 ггг-
0 150000 300000 450000 600000 время, фс
750000 900000 1050000
4
3
2
А
Б
В
Г
Д
Е
Рис. 2. Изменение среднего квадратичного отклонения в ходе молекулярной динамики 48БЯ
35,8 35,6 < 35,4 | 35,2 ! 35
« 34,8 -у
> 34,6
Л 34,4 а.
34,2 34
0
450000 600000 время, фс
900000 1050000
150000
300000
Рис. 3. Изменение радиуса гирации в ходе молекулярной динамики 48БЯ
49000
Н 48000
П е£.. 47000 1 2
д 5 46000
■ Ф
Г * 45000 ; а.
! 44000 х £ 43000 В 42000 41000
щшщщтчпт* г1г"п» г щ 1 жившмшн^ Г711
1ШШдгт1 ншШ г|Ш и л А
ГГ
1
200000
400000
600000 время,фс
800000
1000000
0
Рис. 4. Изменение поверхности, доступной растворителю, в ходе молекулярной динамики 48БЯ
Рис. 5. Изменение среднего квадратичного отклонения в ходе молекулярной динамики для аминокислотных остатков 19, 25, 163, составляющих активный центр 48БЯ
время, фс
Рис. 6. Изменение площади полости активного центра 48БЯ в ходе молекулярной динамики
Рис. 7. Площадь полости активного центра мутанта 48БЯ
является тот факт, что колебания при расчете методами А, Г и Е медленно возрастают на протяжении всего моделирования, хотя для методов А и Г изменение среднего квадратичного отклонения не показывает, что какое-либо подобное изменение начинает происходить. Возрастание среднего квадратичного отклонения в течение молекулярной динамики для метода Е подтверждается медленным увеличением колебаний радиуса гирации. Относительно результатов, полученных методами В и Б, самые лучшие показатели дает метод Д.
Относительно исходной структуры наблюдается расширение поверхности, доступной растворителю (рис. 4). Для методов В (PME), Б расширение поверхности относительно стабильно и колеблется в пределах 42000-44000 Â2, 4400046500 Â2 соответственно. Методы А и Г дают практически идентичные результаты, которые колеблются в пределах 44000-45500 Â2 и по ходу молекулярной динамики имеют тенденцию к уменьшению. Именно эти два метода дают самые лучшие результаты относительно эталонного метода В (PME). Для методов Д и Е расширение поверхности заканчивается примерно на 300 пс моделирования и колеблется в пределах 46000-47000 и 46000-48000 Â2 соответственно. Возрастание колебаний после 900 пс в ходе молекулярной динамики для метода Е согласуется с результатами, полученными при расчете радиуса гирации для данного метода.
В целом наилучшие результаты относительно метода PME показали методы: А - метод с использованием радиуса обрезания 8 Â в программе GPAMM, Г -метод с измененными 1-4 взаимодействиями, как в amber, в программе GPAMM. Эти методы показали лучшие результаты, чем метод расчета с помощью радиуса обрезания потенциала в программе Gromacs-3.3.1 на CPU. Однако по результатам расчета радиуса гирации наилучшим методом расчета относительно PME является метод Д с использованием методики, описанной в [12].
Активный центр
Активный центр представлен тремя аминокислотами (рис. 7). Предполагается, что гистидин принимает непосредственное участие в катализе, глутамин и серин стабилизируют кремнекислородный тетраэдр водородными связями, а аминокислоты, расположенные по бокам трех основных, принимают участие в формировании полости активного центра.
Был произведен расчет среднего квадратичного отклонения для аминокислотных остатков активного центра (рис. 5). Относительно метода В (PME) наилучшие результаты показывают методы Д, Е. Среднее квадратичное отклонение от исходной структуры для методов А, Г медленно возрастает в ходе молекулярной динамики. Был обнаружен интересный факт: для метода Б наблюдаются довольно сильные колебания значений в ходе динамики, а к концу моделирования среднее квадратичное отклонение начинает уменьшаться и приближается к значениям, получаемым методами В, Д, Е. Из двух методов Д и Е наилучшие результаты показывает метод Д, причем среднее квадратичное отклонение для этого метода колеблется гораздо меньше и находится в пределах 0.3-0.6 Â, максимально приближаясь к значениям, полученным методом PME.
Изменение площади треугольника, образованного тремя атомами Nd1, OE1, OG (рис. 7) активного центра, в ходе молекулярной динамики представлено на рис. 6. Относительно эталонного метода В (PME) наилучшие результаты пока-
Табл. 1.
Средние параметры (площадь, расстояние между атомами в аминокислотных остатках), исследованные для каждого из методов
Параметры Методы
А Б В Г Д Е
Площадь, А2 11.139 6.504 7.069 11.202 4.607 5.418
GLN19 -HIS163, А 5.199 9.630 4.572 5.237 4.834 5.499
SER25 - HIS163, А 4.341 4.992 3.675 4.373 3.303 4.033
GLN19 - SER25, А 6.083 5.082 4.088 6.124 2.917 2.858
зывают методы: Б, Д, Е. Для метода Б значения становятся относительно стабильными после 100 пс моделирования и колеблются в пределах 2-15 А2, практически совпадая с данными, получаемыми методом PME. Для методов А, Г площадь медленно возрастает в ходе молекулярной динамики, и эти методы показывают худшие результаты. Самые близкие результаты к данным, получаемым PME, дает метод Е.
В ходе молекулярной динамики был также произведен расчет изменения расстояний между атомами OG 19 глутамина и Nd1 163 гистидина, OE1 25 се-рина и Nd1 163 гистидина, OG 19 глутамина и OE1 25 серина (табл. 1). Наилучшие результаты по всем параметрам табл. 1 относительно эталонного метода В (PME) показывает метод Е.
Для активного центра мутантного белка 4SER, в отличие от данных, полученных нами для всей структуры, наилучшие результаты относительно метода PME показывает метод Е.
Заключение
Было произведено сравнение различных схем расчета электростатических взаимодействий в алгоритмах, написанных для расчета на графических акселераторах. Среди исследованных нами методов относительно метода PME наилучшие результаты при изучении целой структуры мутанта 4SER, растворенного в воде, показали методы: А - метод с использованием радиуса обрезания 8 А в программе GPAMM, Г - метод с измененными 1-4 взаимодействиями, как в amber, в программе GPAMM. При изучении активного центра мутантного белка 4SER наилучшие результаты показал парный метод, основанный на сум-мации Вульфа [13] с радиусом обрезания потенциала 9 А.
Summary
N.I. Akberova, D.A. Alisheva, E.D. Izotova, D.S. Tarasov. Comparison of the Methods of Electrostatic Interactions Treating in Algorithms for Graphical Accelerators.
Molecular dynamics simulations of silicatein alpha cathepsin L chimeras using particle mesh Ewald (PME) and cutoff methods for CPU and pairwise methods for GPU were performed. The results of the pairwise method based on Wolf summation with cutoff 9 А shows best correlation with results of PME method for the active centre of 4SER chimera. For the whole 4SER chimera best results were shown by methods with cutoff 8 А and switching function.
Key words: graphical accelerators, pairwise methods of treating electrostatic interactions, cathepsin L.
Литература
1. Hummer G., Pratt L.R., Garcia A.E. Molecular theories and simulation of ions and polar molecules in water // J. Phys. Chem. A. -1998. -V. 102. - P. 7885-7895.
2. Darden T., York D., Pedersen L. Particle mesh Ewald: an N-log(N) method for Ewald sums in large systems // J. Chem. Phys. - 1993. - V. 98, No 12. - P. 10089-10092.
3. Luty B.A., van Gunsteren W.F. Calculating Electrostatic Interactions using the Particle-Particle Particle-Mesh Method with Nonperiodic Long-Range Interactions // J. Phys. Chem. - 1996. - V. 100. - P. 2581-2587.
4. Steinbach P.J., Brooks B.R. New spherical-cutoff methods for long-range forces in mac-romolecular simulation // J. Comput. Chem. - 1994. - V. 15, No 7. - P. 667-683.
5. Fox T., Kollman P.A. The Application of different solvation and electrostatic models in molecular dynamics simulations of ubiquitin: how well is the X-ray structure "maintained" // Proteins. - 1996. - V. 25, No 3. - P. 315-334.
6. Mark P., Nilsson L. Structure and dynamics of liquid water with different long-range interaction truncation and temperature control methods in molecular dynamics simulations // J. Comput. Chem. - 2002. - V. 23, No 13. -P. 1211-1219.
7. Saito M. Molecular dynamics simulations of proteins in solution: artifacts caused by the cutoff approximation //J. Chem. Phys. - 1994. - V. 101. - P. 4055-4061.
8. Frenkel D., Smit B. Understanding molecular simulation. From Algorithms to Applications. - London: Academic Press, 1996. - 638 p.
9. Leach A.R. Molecular modeling. Principles and applications. - Harlow, England: Pretince-Hall, 2001. - 759 p.
10. ZwartS.F.P., BellemanR.G., Geldof P.M. High-performance direct gravitational N-body simulations on graphics processing units // New Astronomy. - 2007. - No 12. - P. 641650.
11. Yang J., Wang Y., Chen Y. GPU accelerated molecular dynamics simulation of thermal conductivities // J. Comp. Physics. - 2007. - V. 221, No 2. - P. 799-804.
12. Fennell Ch.J., Gezelter J.D. Is the Ewald summation still necessary? Pairwise alternatives to the accepted standard for long-range electrostatics // J. Chem. Phys. - 2006. -V. 124. - P. 234104-1-234104-12.
13. Wolf D., Keblinski P., Phillpot S.R., Eggebrecht J. Exact method for the simulation of coulombic systems by spherically truncated, pairwise 1/r summation // J. Chem. Phys. -1999. - V. 110, No 17. - P. 8254-8282.
14. Zahn D., Schilling B., Kast S.M. Enhancement of the wolf damped coulomb potential: static, dynamic, and dielectric properties of liquid water from molecular simulation // J. Phys. Chem. B. - 2002. - V. 106. - P. 10725-10732.
15. LevittM., HirshbergM., Sharon R., Laidig K.E., Dagett V. Calibration and Testing of a Water Model for Simulation of the Molecular Dynamics of Proteins and Nucleic Acids in Solution // J. Phys. Chem. B. - 1997. - V. 101. - P. 5051-5061.
16. Van der Spoel D., van Maaren P.J., Berendsen H.J.C. A systematic study of water models for molecular simulation: Derivation of water models optimized for use with a reaction field //J. Chem. Phys. - 1998. -V. 108. - P. 10220-10230.
17. Perrera L., Essman U., Berkowitz M.L. Effect of the treatment of long-range forces on the dynamics of ions in aqueous solutions // J. Chem. Phys. - 1995. - V. 102. - P. 450456.
18. Boresch S., Steinhauser O. Rationalizing the effects of modified electrostatic interactions in computer simulations: The dielectric self-consistent field method // J. Chem. Phys. -1999. - V. 111. - P. 8271-8274.
19. Weber W., Hunenberger P.H., McCammon J.A. Molecular Dynamics Simulations of a Polyalanine Octapeptide under Ewald Boundary Conditions: Influence of Artificial Periodicity on Peptide Conformation // J. Phys. Chem. B. - 2000. - V. 104. - P. 36683675.
20. Smith P.E., Pettitt B.M. Peptides in ionic solutions: a comparison of the Ewald and switching function techniques // J. Chem. Phys. - 1991. - V. 95. - P. 8430-8441.
21. York D.M., Darden T.A., Pedersen L.G. The effect of long-range electrostatic interactions in simulations of macromolecular crystals: a comparison of the Ewald and truncated list methods // J. Chem. Phys. - 1993. - V. 99. - P. 8345-8348.
22. De Souza N.O., Orstein R.L. Molecular dynamics simulations of a protein-protein dimer: particle-mesh Ewald electrostatic model yields far superior results to standard cutoff model // J. Biomol. Struct. Dyn. - 1999. - V. 16, No 6. - P. 1205-1218.
23. Тарасов Д.С., Изотова Е.Д., Алишева Д.А. GPAMM [Электронный ресурс]. - 2008. -Режим доступа: www.gpamm.mntech.ru, свободный.
24. FairheadM., Kenneth. A. Johnson, Kowatz T., Stephen A. McMahon, Lester G. Carter, Oke M., Liu H., James H. Naismith, Christopher F. Van der Walle Crystal structure and silica condensing activities of silicatein a-cathepsin L chimeras // Chem. Comm. - 2008 -P. 1765-1767.
25. van der Spoel D., Hess B., LindahlE. Gromacs-3.3.1 [Электронный ресурс]. - Режим доступа: www.gromacs.org.
26. Wang J., Cieplak P., Kollman P.A. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? // J. Comp. Chem. - 2000. - V. 21, No 12. - P. 1049-1074.
27. HyperChem Release 7.0, Computational Chemistry, HyperCube, Inc.: Publication HC70-00-04-00- 2002.
28. Humphrey W., Dalke A., Schulten K. VMD: visual molecular dynamics // J. Mol. Graph. -1996. - V. 14, No 1. - P. 33-38.
29. Pedretti A., Villa L., Vistoli G. VEGA: a versatile program to convert, handle and visualize molecular structure on windows-based PC // J. Mol. Graph. Model. - 2002. - V. 21, No 1. - P. 47-49.
Поступила в редакцию 21.02.08
Акберова Наталья Ивановна - кандидат биологических наук, старший научный сотрудник кафедры биохимии Казанского государственного университета.
E-mail: nakberov@ksu.ru
Алишева Диана Айратовна - студент кафедры биохимии Казанского государственного университета.
E-mail: dalisheva@mntech.ru
Изотова Екатерина Дмитриевна - студент кафедры биохимии Казанского государственного университета.
E-mail: ied@mntech.ru
Тарасов Денис Станиславович - кандидат биологических наук, инженер кафедры биохимии Казанского государственного университета.
E-mail: dtarasov@compnera.com
