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

\

1 7

УДК 612.841.5:532.72:51

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

М.М. Стольниц, А.Н. Башкатов, Э.А. Генина, В.В. Тучин

Саратовский государственный университет, кафедра оптики и биомедицинской физики E-mail: StolnizMM@info.sgu.ru

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

Mathematical Model of Drags and Immersion Liquids Diffusion in Human Ocular Tissues

M.M. Stolnitz, A.N. Bashkatov, E.A. Genina, V.V. Tuchin

Mathematical model of drags and immersion liquids diffusion in human ocular tissues is presented. Analytic expressions are derived for diffusing liquid concentration, which give good approximations for exact solution at small or large time.

Введение

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

В настоящее время существенный прогресс в разработке данных моделей достигнут, например, в области дерматологии [1-4].

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

Другой важной проблемой оптики биотканей является проблема транспорта оптического излучения через фиброзные ткани человеческого тела. Сильное светорассеяние приводит к существенным потерям в данных биотканях зондирующего излучения, снижению разрешающей способности диагностических и эффективности терапевтических методов [5-7]. Управление оптическими па-

раметрами биотканей является решением данной проблемы. Хорошо известно, что рассеяние в биотканях определяется различием показателей преломления структурных компонент биотканей (коллагеновых волокон, митохондрий и т.д.) и окружающей их внутритканевой жидкости. Замещение внутритканевой жидкости иммерсионной жидкостью с показателем преломления большим, чем у внутритканевой жидкости, приводит к существенному снижению светорассеяния. Однако, несмотря на многочисленные исследования, связанные с управлением оптическими параметрами биотканей [5, 8-10], решение проблемы еще далеко от завершения. В настоящее время практически полностью отсутствуют математические модели, описывающие взаимодействие биотканей с осмотически активными иммерсионными жидкостями, и надежные экспериментальные данные по коэффициентам диффузии.

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

1. Математическая модель

Структурно в тканях глаза человека можно выделить три основных слоя [11, 12], характеризующихся различными коэффициентами диффузии. Это конъюнктива (толщина слоя порядка 200 мкм), склера (толщина слоя порядка 500 мкм) и пигментный эпителий (толщина слоя порядка 10 мкм). Однако последний слой намного тоньше второго и может рассматриваться как бесконечно тонкая пленка с эффективным коэффициентом проницаемости.

© М.М. Стольниц, А.Н. Башкатов, Э.А. Генина, В.В. Тучин, 2008

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

и

Рис. 1. Схема слоёв глазной ткани с нанесённым слоем жидкости: 1 - конъюнктива; 2 - склера; 3 - пигментный эпителий. Для наглядности соотношения толщин слоёв не соблюдены

Тогда с точки зрения математической физики задачу можно сформулировать следующим образом. Пусть имеется плоская двухслойная среда, однородная в поперечном направлении, так что задача сводится к одномерной. На верхней границе первого слоя осуществляется обмен со слоем жидкости постоянной концентрации, причем выполняются граничные условия 1-го рода. Систему координат выберем в соответствии с рис. 2. Здесь а,, ^ (у = 1, 2) - координата начала и толщина у-го слоя1 соответственно. Распределение концентрации жидкости и, как функции координаты г и времени г есть решение уравнения диффузии в каждом слое, характеризующемся коэффициентом диффузии В,. На границах между слоями выполняются условия сопряжения, выражающие равенство концентраций и потоков. На нижней границе второго слоя выполняются граничные условия 3-го рода - происходит обмен со средой с нулевой концентрацией диффундирующей жидкости (физически это соответствует постоянному оттоку жидкости от границы и, следовательно, отсутствию её обратного потока в биоткань).

1 Практически всегда можно положить а1 = 0, как

это и изображено на рисунке, однако для упрощения это

сделано в окончательных выражениях.

—— и В 1 —►

а3 и В2 К32 2

и = 0

Рис. 2. Система координат и обозначения для краевой задачи

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

и(г,г) = и. (г,г), при а. < г < ау+1.

Тогда краевая задача для нее включает: уравнение диффузии

ди, (г, г) 2 д2и, (г,г) 2

=«5-Й-2 • 5 = В • I=и ; ш

граничные условия

и1 (а1’ г) = и0’

г 2 ди2 (г,г) , ч

82 2д( ) + К 32 и 2 (г, г)

V дг У

= 0 (2)

г=аз

где и0 - постоянная концентрация жидкости за границами исследуемой среды; к32 - коэффициент проницаемости на границе 2-го слоя и окружающей среды, который определяется через коэффициент диффузии и толщину третьего слоя: к32 = 532/d3;

условия сопряжения на границе слоев

и2 (а2,г) = и1 (а2,г), ди1 (г, г)

5

дг

= 52

дг

(3)

г=а2

г=а2

начальные условия

и ( г, г ) |г=0 = 0 при а1 < г < а3 . Введём безразмерные переменные:

г - а.

(4)

гу =

' 542

. d.

V , У

иу (г,.г,)

- при а, < г < а,+,,

и

У = 1,2

(5)

^ = а,+1- а,

2

При этом краевая задача (1) - (4) запишется в виде

дuj (^, tj)=дЧЛ

u1 (0, ^) = 1,

2 ди2 (22, t2 ) ? 1 \

52 - + K5з2u2 ( ^2, t2 )

аг2

1 = 1,2, (1')

(2')

= 0,

22 =1

ди2 ( 22 , t2 )

= од

ди1 (^ t1)

^2 22=° ^

иі ( )| t. = ° = 0’

(з')

(4')

где

К«32 = К 32 4 = «^23. t! = МгА,

«2

ОД. = 5,2 Д., О. = «-, 4 ^

1 1 . 1 О ’ 1 Д

1 1

(6)

2. Результаты и их обсуждение

Аналитическое решение краевой задачи (1') - (4') получим с помощью преобразования Лапласа [13, 14]. При начальных условиях вида (4') применение преобразования Лапласа приводит к общему выражению, из которого можно получить виды решения, хорошо сходящиеся или на больших, или на малых временах. Решения первого вида получаются с помощью теории вычетов и представляют собой ряды Фурье, а с помощью разложения функции-образа в ряд и обратного преобразования Лапласа его членов получаются решения в виде кратных рядов экспонент и функций ошибок. Хотя такой ряд громоздок и мало пригоден для практических расчетов, уже первые его члены дают достаточную точность как раз в том интервале значений г, где плохо сходятся ряды Фурье. Таким образом, комбинация решений обоих видов охватывает весь интервал времени проведения эксперимента.

Применим преобразование Лапласа

tdt

к соотношениям (1')-(4'). Пользуясь правилами работы с Ь-образами [15], сведем краевую задачу к системе обыкновенных дифференциальных уравнений:

Д) 2 / \ Д, і—

----«1и,ь(2И- «1 =57^ 1=1-2

и. (°)= -

и2Ь ( 0) = и1Ь (1) , Ди2Ь ( г2 )

(7)

Ди1Ь ( г1 )

22 =0

Ді1

21=1

!- \ , 1 Ди2Ь (22 ) 1

и2Ь ( 22 ) + о ,

КО32

= 0

22 = 1

Начальное условие при этом автоматически учтено.

Будем искать решение (7) в виде

и1Ь (^ («і2і) + ві сЬ («і2і). (8)

Подставляя (8) в (7), после ряда преобразований получаем:

1 «12 ( «2 ) СЬ ( «1 )-^2 ( «2 ) ЄІ1 ( «21 )

,( ) =

,( 22 ) =

5 К«з2 Я1 ( «1, «2 ) + «2 ^2 ( «1, «2 )

«2 ^ ( «22 ) + «2 сЬ ( «22 )

5 К«32^1 ( «1, «2 ) + «2К2 («1, «2 )’

(9)

где

^ (42) = к532 (<?2) + 42сЬ (),

^2 ( 42 ) = к532 сЬ (42 ) + 42 ^ (42 ) .

4г1 = 41 (г1 -1),

Щ (41.42 ) = 512сЬ (41) (42) + вЬ (41) А (42),

к2 (41,42 ) = 512сЬ (41) сЬ (42) + (41)(42),

4г2 = 42 (1- г2).

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

I \ - к532

и1 (*1,г1 ) =1 —„ ^ ^ , г +

+21

К«32 ( ОД21 + 1) + ОД

Р1 (^1т , 21 ) ЄХР (^\ )

(10)

^1тб (^1т )

и

т=1

0

р1 (Мш. Z) = [ к5з2 COS ( Si2 d 2iM-1m ) +

+812d21^1m Sin (812d21^1m )] Sin |>1m (1 - Z )] -

-812 [ k832 Sin ( 812d21^1m ) +

+ 812d21^1m COS (812d21^1m )] COS [Mb» (1 “ Z )_,

P2 ( M1m • Z2 ) = k832 Sin [812d21^1m (1 “ Z2 )] +

+ §12d21^1m COS [812d21^1m (1 “ Z,

Q ( M-1m ) = Q1 (ц1т ) — 812d21M1mQ2, 1m )

Q1 (^1m ) = [ K832 C1 + 8d21 ) + 8d21 ] X

X COS ( Ц 1m ) COS ( 812 d 21Ц 1m )“

“ 812 [ K832 (1 + d21 ) + d21 ] X X Sin (Ц 1m ) Sin (812d, М 1m )

Q2 (Mlm ) = 812 (1 + d21 ) Sin (^1m ) COS (812d21M1m ) + + (1 + 8d21 ) COS ( M"1m ) Sin ( 812 d21^1m ,

где M1m, m = 1,2,3,... - положительные корни уравнения

k832 [812 cos (ц1 ) Sin (812 d 21ц1 ) +

+ Sin (Ц1) cos (812d 21Ц)] + (12)

+ 812 d21Ц1 [ 812 COS (ц1 ) COS ( 812 d21 М-1 ) “

- sin (М1) sin (812d21М1)] = 0.

При к832 =~ решение (10) совпадает с решением из [14, гл. XII, §8]. Первые члены в формулах (10), очевидно, представляют стационарное решение, к которому приближаются распределения концентрации жидкости в обоих слоях при t ^ ^. Следовательно, с течением времени вклад каждого ряда в значения концентраций уменьшается, и для достижения необходимой точности приближения требуется все меньше первых членов ряда. Уже при t1 > 1 для практических расчетов достаточно первых пяти членов. С другой стороны, при t1 » 0 даже для получения общего вида распределений концентрации жидкости, качественно совпадающего с реальным, требуется больше 1000 членов ряда. Результаты расчетов по формулам (10) приведены на рис. 3. Следует отметить не имеющие физического смысла осцилляции решения при t = 0 и недостаточно большом числе членов ряда Фурье.

Рис.3. Распределение концентрации диффундирующей жидкости при различных значениях г, рассчитанное по формулам (10) (учитываются первые 10 членов ряда Фурье): 1 - г\ = 0; 2 - г\ = 0.1; 3 - г\ = 1; 4 - г\ = да

Чтобы получить компактное решение, дающее хорошую точность на малых временах, нужно представить (9) в виде ряда экспонент. Ограничиваясь членами первого порядка малости, имеем:

I \ 1 J 812 “1

U1L (Z Н-

exp I-^ (2 - z1 )Vj[ +

(

+ exp I -

d „ d —z1 + 2v 81 1 8,

+ exp

. 81

vj J

(d d Л

f (2 - z ) + 2 d2

81 82 у

d1

+

vj!

+

+eXP |- Z1^J J

2k81

(13)

32

X

k832 + Vj

exp 810 -1

X

8-(2-Z1)+28. 4-j j+

f

+

812 +1

eXP I-

dL + 2^

81 1 8.

2 у d Л

'2 У

vj!

.( 2 2 )

28,

( 812 + 1) 5 + ехр\ -

ехр<

Гй1 й2 Л

— + — 22

81 8 2 У

+

г й1 й 2 Л

— + 22

81 82 2 у

л/5!

+ ехр | -2 к8'

л/л-

32

К8 32 + л/5

К8 32 = К8 32 8~ й 2

ехр

Т1 + ( 2 - 22 )

81 8 2

л/5!

Использование таблиц обратного преобразования Лапласа (см. [15, приложение, формулы (168), (170)) дает:

( \ 812 -1 Щ (21, г1 )-

812 +1

еГс

2 - 21

- еГС

1

г \_

21 + 2812й21

V

У

+

/ \ /

егЪ 21 - еТс

26 У V

2 - 21 + 2812й21

+

+2ехр

х<!ехр

2

к832

V 812й21 у

^1 + 2К832

х

К832 ( 2-21) 812й21

еГС

2 21 +2812й21 + К832

812 -1 +——ехр

812 +1

к83221 812й21

812й21

+

еГС

>( 22, Ь1)-

_281 812 +1

ег&

< ег^с

21 + 2812й21 + к832 21ь1 812й21

Л

(14)

2л/г1 ,

1 + 812 й21 (2 - 22)

+ 2ехр

2

К832Л V 812 й21 у

2д/г1

к8.0

812 й21

+

I +

[1 + 812 й21 ( 2 - 22 ) ]

х

х eгfc

1 + 812 й21 ( 2 22 )+ к832 ^1

2у[г1

812 й21

Анализ выражений (14) показывает, что они дают точное решение при г1 = 0. Действительно, в этом случае все члены обеих формул, кроме одного, равны 0 при всех значениях 21, 22. Что же касается третьего члена в и1 (21, г1), то при 21 = 0 в его аргументе появляется неопределенность «0/0», и его значение зависит от порядка перехода к пределам 21 ^ +0, г1 ^ +0. А это полностью совпадает с начальными условиями: нулевое значение концентрации на всем протяжении исследуемой ткани и единичное - на границе первого слоя. Если же ь1 > 0, но всё еще мало (г1 << шт{2,812й21}), то щ (21, г1) - 0, при

21 >> Ь1 и и1 (21, Ь1) - 1 в противоположном случае. Таким образом, наблюдается резкий «диффузионный фронт», который становится более пологим при больших г1). Результаты расчетов по формулам (14) приведены на рис. 4.

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

В качестве примера рассмотрим три лекарственных препарата: инулин, надолол и тинолол. В таблице приведены значения ко-

Рис. 4. Распределение концентрации диффундирующей жидкости при различных значениях г, рассчитанное по формулам (14): 1 - г1 = 0; 2 - г1 = 0.01; 3 - г1 = 0.09; 4 - г1 = 0.09 (кривая 4 рассчитана по формуле (10))

и

г

\

\

\

\

и

Г

\

Диффузионные характеристики тканей глаза для трёх лекарственных препаратов

Г лазная ткань Диффун- дирующее вещество Коэффициент проницаемости P, м/с Коэффициент распределения K Коэффициент диффузии D, м2/с

Конъюк- тива Инулин 3.8Х10-8 1 7.8Х10-12

Надолол 5.3Х10-7 0.23 4.61Х10-10

Т инолол 5.2Х10-7 1.61 6.46Х10-11

Склера Инулин 9.0Х10-8 1 4.5Х10-11

Надолол 3.9Х10-7 0.23 8.48Х10-10

Т инолол 4.1Х10-85 1.61 1.27Х10-10

эффициентов проницаемости и «коэффициентов распределения» конъюнктивы и склеры согласно обзору [16]. В свою очередь, коэффициент проницаемости P связан с коэффициентом диффузии D выражением Р = = KD/ l [2, 17], где К — коэффициент распределения и l — толщина слоя. Рассчитанные коэффициенты диффузии данных лекарственных препаратов приведены в последнем столбце таблицы. При этом безразмерный параметр к532 может принимать значения из диапазона 20-50. Это позволяет рассчитать, например, время выхода концентрации диффундирующего вещества на стационар £стац, определяемое как промежуток времени от начала процесса диффузии, за которое распределение концентрации становится близким к стационарному с некоторой заданной точностью, например, неотличимо визуально на графике. Для безразмерной величины ?1стац расчёты дают значение «3.5, что соответствует времени £стац « 1.8 х 104, 300 и 2150 с для инулина, надолола и тинолола соответственно.

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

Работа выполнена при финансовой поддержке РФФИ (грант № 06-02-16740-а) и U.S. Civilian Research and Development Foundation for the Independ-

ent States of the Former Soviet Union (грант № REC-006/SA-006-00).

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

1. Томашевич М.С., Воротникова Т.В., Быков В.А. Моделирование систем доставки биологически активных веществ в кожу и подлежащие ткани // Биомедицинские технологии и радиоэлектроника. 2006. №3. С.31-37.

2. Schaefer H., Redelmeier T.E. Skin Barrier: Principles of percutaneous absorption. Basel: Karger, 1996. 322 p.

3. Watkinson A.C., Brain K.R. Basic mathematical principles in skin permeation // J. Toxicology: Cutaneous and Ocular Toxicology. 2002. Vol.21, №4. P.371—402.

4. Yamashita F., Bando H., Koyama Y. et al. In vivo and in vitro analysis of skin penetration enhancement based on a two-layer diffusion model with polar and nonpolar routes in the stratum corneum // Pharm. Res. 1994. Vol.11, №2. P.185—191.

5. Bashkatov A.N., Genina E.A., Tuchin V.V. Optical immersion as a tool for tissue scattering properties control // Perspectives in Engineering Optics / Eds. K. Singh, V.K. Rastogi. New Delhi, India: Anita Publications, 2002. P.313—334.

6. Tuchin V.V. Optical clearing in tissues and blood. Bellingham, WA: SPIE Press, 2005. 256 p.

7. Tuchin V.V. Tissue optics: Light scattering methods and instruments for medical diagnosis. Bellingham, WA: SPIE Press, 2007. 840 p.

8. Башкатов А.Н., Генина Э.А., Синичкин Ю.П., Кочубей В.И., Лакодина Н.А., Тучин В.В. Определение коэффициента диффузии глюкозы в склере глаза человека // Биофизика. 2003. Т. 48, №2. С.309—313.

9. Bashkatov A.N., Genina E.A., Sinichkin Yu.P., Kochubey V.I., Lakodina N.A., Tuchin V.V. Glucose and mannitol diffusion in human dura mater // Biophys. J. 2003. Vol.85, №5. P.3310—3318.

10. Tuchin V.V. Optical clearing of tissues and blood using the immersion method // J. Phys. D: Appl. Physics. 2005. Vol.38, №15. P.2497—2518.

11. Hammer M., Roggan A., Schweitzer D., Muller G. Optical properties of ocular fundus tissues - an in vitro study using the double-integrating-sphere technique and inverse Monte Carlo simulation // Phys. Med. Biol. 1995. Vol.40, №6. P.963—978.

12. Nemati B., Dunn A., Welch A.J., Rylander III H.G. Optical model for light distribution during transscleral cyclophotoco-agulation // Appl. Opt. 1998. Vol.37, № 4. P.764—771.

13. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1972. 736 с.

14. Карслоу Г., Егер Д. Теплопроводность твердых тел / Пер. с англ. М.: Наука, 1964. 488 с.

15. Дёч Г. Руководство к практическому применению преобразования Лапласа и Z-преобразования / Пер. с англ. М.: Наука, 1971. 288 с.

16. Prausnitz M.R., Noonan J.S. Permeability of cornea, sclera, and conjunctiva: a literature analysis for drug delivery to the eye // J. Pharm. Sci. 1998. Vol.87, №12. P.1479—1488.

17. Котык А., Яначек К. Мембранный транспорт. Междисциплинарный подход / Пер. с англ. М.: Мир, 1980. 342 с