Line Tension of Pore Edge in Membrane on Solid Support

Full Text

ВВЕДЕНИЕ

Бислойные липидные мембраны являются чрезвычайно слабопроницаемыми тонкими пленками [1]. Высокие барьерные свойства мембран используются живыми клетками для поддержания постоянного состава цитоплазмы, отличного от состава окружающей среды, и состава клеточных органелл, отличного от состава цитоплазмы. Малая толщина мембран (~ 4 нм) позволяет поддерживать значительные градиенты концентраций различных веществ и электрических потенциалов [2]. В норме транспорт необходимых клетке веществ через мембрану осуществляется посредством специализированных каналов или насосов [3, 4] либо посредством топологической перестройки мембраны в процессах экзо- и эндоцитоза [5, 6], а также фагоцитоза [7]. Резкое неконтролируемое повышение проницаемости мембраны вследствие образования в ней сквозной поры, как правило, приводит к гибели клетки [8]. Однако понимание механизмов контролируемого образования пор в липидных мембранах является ключом к разработке новых систем трансфекции и адресной доставки лекарств, технологии создания гибридных клеток и т.д.

Впервые образование пор в тонких бесструктурных пленках было теоретически описано в работе Дерягина и Гутопа [9]. В данной теории предполагается, что энергия пленки с порой радиуса R может быть записана в виде: E(R) = 2pRg – pR²s. Первое слагаемое, пропорциональное периметру поры (2pR) с коэффициентом пропорциональности g, называемым линейным натяжением [9], соответствует изменению внутренней энергии пленки при образовании поры. Это слагаемое учитывает различие энергии материала пленки на кромке поры и вдали от нее. Второе слагаемое, пропорциональное площади поры (pR²) с коэффициентом пропорциональности s, называемым латеральным натяжением, учитывает работу внешних сил. Соответственно, латеральное натяжение является для мембраны внешним параметром, определяющимся конкретной экспериментальной системой; как правило, латеральное натяжение точно контролируется в эксперименте [10–12]. Линейное натяжение, напротив, является параметром, который определяется физико-химическими свойствами материала пленки, а также детальной структурой кромки поры. Линейное натяжение не может контролироваться в эксперименте, однако может быть измерено [13–15]. В модели Дерягина–Гутопа энергия, как функция радиуса, имеет максимум при так называемом критическом радиусе R* = g/s. При 0 < R < R* энергия понижается при уменьшении радиуса; при R > R* энергия бесконечно понижается при увеличении радиуса. Таким образом, состояние пленки без поры (с порой нулевого радиуса) оказывается метастабильным, а основному состоянию системы соответствует пора бесконечного радиуса, т.е. полное отсутствие пленки. Величина энергии DE = E(R*) = pg²/s определяет высоту энергетического барьера, который необходимо преодолеть для перехода из метастабильного (R = 0) в основное (R → ∞) состояние. Чем больше латеральное натяжение и чем меньше линейное натяжение, тем меньше высота энергетического барьера DE.

При R = R* энергия максимальна, и пленка с порой находится в состоянии равновесия, хотя и неустойчивого. Соответственно, в самом положении равновесия потенциальная сила F = –dE(R*)/dR равна нулю. В окрестности критического радиуса эта сила мала, и изменения радиуса под действием потенциальной силы F должны быть относительно медленными. Такое поведение системы регистрировалось экспериментально методом атомной силовой микроскопии (АСМ) [16]. Толстая полимерная пленка, к которой было приложено некоторое латеральное натяжение, наносилась на твердую подложку. В пленке лазерным лучом прожигались отверстия различного радиуса. Отверстия малого радиуса спонтанно закрывались, отверстия большого радиуса увеличивались. При некотором промежуточном радиусе (~ 1 мкм), соответствующем критическому радиусу данной системы, размер отверстия оставался практически постоянным в течение длительного времени (~ 15 мин) [16]. Таким образом, при известном латеральном натяжении s, путем измерения критического радиуса системы R* = g/s, данный метод может использоваться для определения величины линейного натяжения кромки поры g. Для модельных липидных мембран при типичном латеральном натяжении s ~ 1 мН/м [17, 18] и линейном натяжении g ~ 15 пН [13] величина критического радиуса составляет R* ~ 15 нм. Соответственно, такие отверстия не могут прожигаться оптическим лазером. Однако поры такого радиуса, в принципе, могут быть сформированы в липидном бислое, нанесенном на твердую подложку, либо путем нормального продавливания мембраны зондом АСМ подходящего размера в режиме измерений «сила – расстояние» (force – distance) [19, 20], либо путем латерального сдвига липида малым зондом [21].

Состояние мембраны в экспериментах АСМ принципиально отличается от свободной мембраны, например, гигантской однослойной везикулы, наличием твердой слюдяной подложки. Подложка автоматически делает мембрану асимметричной, поскольку с подложкой взаимодействует только один липидный монослой из двух. Ранее было показано [22], что такая асимметрия может изменять среднюю конфигурацию границы жидко-упорядоченных липидных доменов с окружающей жидко-неупорядоченной мембраной, вплоть до практически полного исчезновения одной из двух возможных равновесных конфигураций границы, имеющихся в свободной симметричной мембране. Кроме того, теоретически было предсказано влияние твердой подложки на адсорбцию на мембрану амфипатических пептидов и их латеральное взаимодействие [23]. Возможное влияние подложки на линейное натяжение кромки поры, насколько нам известно, ранее не исследовалось. В настоящей работе проведен расчет линейного натяжения кромки поры в мембране на подложке для различных режимов взаимодействия мембраны с подложкой. Показано, что влияние подложки достаточно мало: линейное натяжение поры в мембране с бесконечно сильной адгезией к подложке отличается от линейного натяжения поры в свободной мембране приблизительно на 3.5%.

ПОСТАНОВКА И РЕШЕНИЕ ЗАДАЧИ

Кромка сквозной поры в липидном бислое характеризуется сильными деформациями липидного материала [24]. Для расчета энергии деформаций мы используем теорию упругости липидных мембран, исходно разработанную Хаммом и Козловым [25], позднее обобщенную с учетом дополнительных деформационных мод [26]. В рамках этой теории упругости деформации считаются малыми, и энергия вычисляется в квадратичном порядке по ним. Средняя ориентация липидных молекул описывается векторным полем единичных векторов n, называемых директорами. Поле директоров задается на поверхности, называемой нейтральной, расположенной внутри липидного монослоя на границе полярных головок и гидрофобных хвостов липидов [27]. Форма нейтральной поверхности задается векторным полем единичных нормалей N к ней. Будем считать, что твердая подложка и мембрана на ней горизонтальны. Ближайший к подложке липидный монослой будем называть нижним, и параметры, относящиеся к нему, будем обозначать индексом «l». Удаленный от подложки липидный монослой будем называть верхним, и относящиеся к нему параметры будем обозначать индексом «u». Введем декартову систему координат Oxyz с началом координат O, расположенным на поверхности подложки, и осью Oz, направленной перпендикулярно подложке в сторону мембраны (рис. 1). В этой системе координат формы нейтральных поверхностей верхнего и нижнего монослоев, а также межмонослойной поверхности могут быть заданы функциями z-координат их точек: H_u(x, y), H_l(x, y), M(x, y) соответственно. Поверхностная плотность энергии деформированного монослоя может быть записана в виде [26]:

$\begin{array}{l}w=\frac{k_c}{2}{\left[div\left(\mathit{n}\right)+J_0\right]}^2-\frac{k_c}{2}{J}_0^2+\frac{k_t}{2}{\mathit{t}}^2+\frac{k_a}{2}{\alpha }^2+\\ +\frac{\sigma }{2}{\left[\mathit{g}\mathit{r}\mathit{a}\mathit{d}\left(H\right)\right]}^2+k_{G}K+\frac{k_{rot}}{2}{\left[\mathit{r}\mathit{o}\mathit{t}\left(\mathit{n}\right)\right]}^2,\end{array}$ (1)

где k_c, k_t, k_a, k_G, k_rot – модули поперечного изгиба, наклона, латерального растяжения, гауссовой кривизны и кручения соответственно; J₀ – спонтанная кривизна монослоя; t = n – N – вектор наклона; a – относительное изменение площади нейтральной поверхности; $K=\frac{\partial n_x}{\partial x}\frac{\partial n_y}{\partial y}-\frac{\partial n_x}{\partial y}\frac{\partial n_y}{\partial x}$ – гауссова кривизна (n_x, n_y – соответствующие проекции директора); H(x, y) – z-координаты точек нейтральной поверхности монослоя; s – латеральное натяжение. Первое слагаемое соответствует деформации поперечного изгиба, который для типичных липидных мембран является самой мягкой деформационной модой [28, 29]. При отсутствии деформации наклона, т.е. при совпадении директора и нормали, n = N, div(n) = div(N) = –J, где J – геометрическая кривизна поверхности, рассчитываемая как сумма принципиальных кривизн C₁, C₂ поверхности: J = C₁ + C₂. Для сквозной поры, форма кромки которой может приближенно рассматриваться как половина тороидальной поверхности, C₁ ~ 1/h – меридиональная кривизна, C₂ ~ 1/r – экваториальная кривизна, где h ≈ 1.5 нм – толщина гидрофобной части монослоя (расстояние между межмонослойной и нейтральной поверхностями) (рис. 1), r – радиус поры в экваториальной плоскости. В окрестности критического радиуса (~ 15 нм), r >> h, следовательно, C₁ >> C₂, и экваториальной кривизной кромки поры можно пренебречь. В этом случае система становится эффективно одномерной: все величины зависят лишь от координаты, перпендикулярной линии кромки поры, а вдоль линии кромки система обладает трансляционной симметрией. Направим ось Ox декартовой системы координат перпендикулярно линии кромки поры, а ось Oy – параллельно ей (рис. 1). Во всех дальнейших расчетах кромка поры предполагается одномерной, т.е. обладающей трансляционной симметрией вдоль оси Oy. В такой трансляционно симметричной системе, описываемой в рамках линейной теории упругости, rot(n) = 0, $K=\frac{\partial n_x}{\partial x}\frac{\partial n_y}{\partial y}-\frac{\partial n_x}{\partial y}\frac{\partial n_y}{\partial x}=\frac{\partial n_x}{\partial x}\cdot 0-0\cdot 0=0$, т.е. деформации кручения и гауссовой кривизны не дают вклада в упругую энергию липидного монослоя. Кроме того, векторные величины могут быть заменены их проекциями на ось Ox: n → n_x = n; N → N_x = N, t = n – N → t = n – N. Дополнительно, div(n) → dn/dx; grad(H) → dH/dx. Таким образом, поверхностная плотность энергии деформации липидного монослоя может быть записана в следующем виде:

$\begin{array}{l}w=\frac{k_c}{2}{\left(n^{\text{'}}+J_0\right)}^2-\frac{k_c}{2}{J}_0^2+\\ +\frac{k_t}{2}{\left(n-N\right)}^2+\frac{k_a}{2}{\alpha }^2+\frac{\sigma }{2}{\left(H^{\text{'}}\right)}^2.\end{array}$ (2)

Здесь и всюду ниже штрих обозначает производную по координате x.

 

Рис. 1. Схематическое изображение одномерной кромки поры в мембране на подложке. Подложка показана наклонной штриховкой. Ось Oz декартовой системы координат перпендикулярна плоскости подложки; ось Ox перпендикулярна линии кромки поры. Форма нейтральных поверхностей верхнего и нижнего монослоев, а также межмонослойной поверхности описывается функциями H_u(x), H_l(x), M(x) соответственно. Толстой красной линией показана нейтральная поверхность вертикального монослойного участка; ее форма описывается функцией H_v(z). Нейтральная поверхность вертикального монослойного участка непрерывно сопрягается с нейтральной поверхностью нижнего монослоя в точке {X₁, Z₁}, верхнего монослоя – в точке {X₂, Z₂}. Толщина гидрофобной части монослоя h.

 

Модуль всестороннего сжатия липидных мембран очень велик, ~ 10¹⁰ Дж/м³ [30]. Можно считать, что гидрофобная часть липидного монослоя локально объемно несжимаема, т.е. что объем элемента монослоя не изменяется при деформации. С требуемой точностью условие локальной объемной несжимаемости для верхнего и нижнего монослоев имеет вид [25, 26]:

$\begin{array}{l}{H}_u-M=h-\frac{h^2}{2}{n^{\text{'}}}_u-h{\alpha }_u,\\ M-H_l=h-\frac{h^2}{2}{n^{\text{'}}}_l-h{\alpha }_l.\end{array}$ (3)

Из этих условий возможно выразить a_u, a_l через n_u, n_l, H_u, H_l, M. Кроме того, с требуемой точностью N_u ≈ H_u′, N_l ≈ –H_l′. Окончательно энергия деформации мембраны записывается в виде:

$\begin{array}{l}W=\int dx\left\{\frac{k_c}{2}{\left({n^{\text{'}}}_u+J_0\right)}^2-\frac{k_c}{2}{J}_0^2\right.+\frac{k_t}{2}{\left(n_u-{H^{\text{'}}}_u\right)}^2+\\ \left.+\frac{k_a}{2h^2}{\left(h-\frac{h^2}{2}{n^{\text{'}}}_u+M-H_u\right)}^2+\frac{\sigma }{2}{\left({H^{\text{'}}}_u\right)}^2\right\}+\\ +\int dx\left\{\frac{k_c}{2}{\left({n^{\text{'}}}_l+J_0\right)}^2-\frac{k_c}{2}{J}_0^2\right.+\frac{k_t}{2}{\left(n_l+{H^{\text{'}}}_l\right)}^2+\\ \left.+\frac{k_a}{2h^2}{\left(h-\frac{h^2}{2}{n^{\text{'}}}_l-M+H_l\right)}^2+\frac{\sigma }{2}{\left({H^{\text{'}}}_l\right)}^2\right\}.\end{array}$ (4)

Интегрирование ведется по координате x; таким образом, выражение (4) дает энергию, отнесенную к единице длины вдоль оси Oy, т.е. вдоль кромки поры.

Необходимо учесть, что нижний монослой горизонтального бислоя взаимодействует с твердой подложкой. В работе [22] рассматривалось влияние подложки на конфигурации границ упорядоченных доменов в мембране. В ней учитывались электростатические, ван-дер-ваальсовы и гидратационные взаимодействия. Полная сила взаимодействия мембраны с плоской заряженной подложкой, отнесенная к единице площади, записывалась в виде [31, 32]:

$f=\frac{2{\sigma }_m{\sigma }_p}{{\epsilon }_0{\epsilon }_w}{e}^{-\frac{z}{l_D}}-\frac{H_a}{12\pi l_D{z}^2}+P_0{e}^{-\frac{z}{l_h}}$, (5)

где s_m – поверхностная плотность электрического заряда мембраны; s_p – поверхностная плотность электрического заряда подложки; e₀ – диэлектрическая проницаемость вакуума; e_w – диэлектрическая проницаемость воды; l_D ≈ 1 нм – длина Дебая в 100 мМ бинарном электролите; H_a ≈ (3 ÷ 10) × 10^–21 Дж – константа Гамакера; P₀ ≈ 4 × 10⁷ ÷ 10¹⁰ Па – расклинивающее давление; l_h ≈ 0.1 – 0.3 нм – характерная длина гидратационного отталкивания [33]. Первое слагаемое в выражении (5) соответствует электростатическим взаимодействиям, второе слагаемое – ван-дер-ваальсовым, третье слагаемое – гидратационным. Энергия взаимодействия мембраны и подложки, отнесенная к единице площади, может быть получена путем интегрирования силы (5) по координате z:

$\begin{array}{c}{w}_s\left(H_l\right)=\underset{\infty }{\overset{H_l-z_0}{\int }}fdz=\frac{2l_D{\sigma }_m{\sigma }_p}{{\epsilon }_0{\epsilon }_w}{e}^{-\frac{H_l-z_0}{l_D}}-\\ -\frac{H_a}{12\pi l_D\left(H_l-z_0\right)}+P_0{l}_h{e}^{-\frac{H_l-z_0}{l_h}}\end{array}$, (6)

где (H_l – z₀) – расстояние от подложки до полярных групп липидов нижнего монослоя мембраны; z₀ ≈ 0.7 нм – толщина полярной части липидного монослоя [27]. При физически разумных значениях параметров эта энергия имеет минимум. Положение минимума H₀ не может быть найдено аналитически, поэтому мы находили его численно при заданных численных значениях параметров. Разлагая энергию в окрестности точки H₀ до второго порядка, найдем эффективный квадратичный потенциал, в котором находится нейтральная поверхность нижнего монослоя мембраны:

$w_s\approx \frac{K_s}{2}{\left(H_l-H_0\right)}^2$. (7)

При типичных значениях параметров s_m = 0, s_p = –0.032 Кл/м² [31], e_w = 81, l_D = 1 нм, l_h = 0.2 нм, P₀ = 4 × 10⁷ Н/м² [33], H_a = 5 × 10^–21 Дж [31], для крутизны потенциала получается значение K_s = 0.3 × 10¹⁵ Дж/м⁴ = 0.08 k_BT/нм⁴ и H₀ ≈ 1.7 нм. Потенциал вида (7) налагает граничное условие на форму нейтральной поверхности нижнего монослоя: необходимо, чтобы H_l(x → ∞) → H₀. При этом конкретное значение величины H₀ не имеет особого физического смысла; H₀ может быть сделано равным нулю путем параллельного переноса системы координат вдоль оси Oz. В дальнейшем для простоты будем считать, что H₀ = 0 и w_s = K_s(H_l)²/2. Значение K_s = 0 соответствует предельному случаю свободной мембраны, K_s = ∞ – мембране с бесконечно сильной адгезией к подложке. В последнем предельном случае H_l ≡ 0. С учетом взаимодействия нижнего монослоя горизонтального бислойного участка с подложкой, функционал энергии горизонтального бислойного участка записывается в виде:

$\begin{array}{c}{W}_h=\underset{X_2}{\overset{\infty }{\int }}dx\left\{\frac{k_c}{2}{\left({n^{\text{'}}}_u+J_0\right)}^2-\frac{k_c}{2}{J}_0^2\right.+\frac{k_t}{2}{\left(n_u-{H^{\text{'}}}_u\right)}^2+\left.\frac{k_a}{2h^2}{\left(h-\frac{h^2}{2}{n^{\text{'}}}_u+M-H_u\right)}^2+\frac{\sigma }{2}{\left({H^{\text{'}}}_u\right)}^2\right\}+\\ +\underset{X_1}{\overset{\infty }{\int }}dx\left\{\frac{k_c}{2}{\left({n^{\text{'}}}_l+J_0\right)}^2-\frac{k_c}{2}{J}_0^2\right.+\frac{k_t}{2}{\left(n_l+{H^{\text{'}}}_l\right)}^2+\left.\frac{k_a}{2h^2}{\left(h-\frac{h^2}{2}{n^{\text{'}}}_l-M+H_l\right)}^2+\frac{\sigma }{2}{\left({H^{\text{'}}}_l\right)}^2+\frac{K_s}{2}{H}_l^2\right\},\end{array}$ (8)

где интегрирование ведется по поверхностям верхнего и нижнего горизонтальных монослоев соответственно. Вариация этого функционала по функциям n_u, n_l, H_u, H_l, M приводит к пяти уравнениям Эйлера–Лагранжа.

$\begin{array}{l}\left(l^2+\frac{A{h}^2}{4}\right){n^{\text{'}\text{'}}}_u-n_u+\left(1+\frac{A}{2}\right){H^{\text{'}}}_u-\frac{A}{2}{M}^{\text{'}}=\mathrm{0,}\\ \left(l^2+\frac{A{h}^2}{4}\right){n^{\text{'}\text{'}}}_l-n_l-\left(1+\frac{A}{2}\right){H^{\text{'}}}_l+\frac{A}{2}{M}^{\text{'}}=\mathrm{0,}\\ \left(1+\frac{A}{2}\right){n^{\text{'}}}_u-\left(1+s\right){H^{\text{'}\text{'}}}_u+\frac{A}{h^2}\left(H_u-M-h\right)=\mathrm{0,}\\ \left(1+\frac{A}{2}\right){n^{\text{'}}}_l+\left(1+s\right){H^{\text{'}\text{'}}}_l+\frac{A}{h^2}\left(M-H_l-h\right)-k_s{H}_l=\mathrm{0,}\\ \left({n^{\text{'}}}_u-{n^{\text{'}}}_l\right)+\frac{2}{h^2}\left(H_u+H_l-2M\right)=\mathrm{0,}\end{array}$ (9)

где A = k_a/k_t, k_s = K_s/k_t. Эта система уравнений линейная с постоянными коэффициентами. Из последнего уравнения системы (9) выразим форму межмонослойной поверхности:

$M=\frac{h^2}{4}\left({n^{\text{'}}}_u-{n^{\text{'}}}_l\right)+\frac{1}{2}\left(H_u+H_l\right)$. (10)

Далее введем новые функции a(x) = n_u(x) + n_l(x); b(x) = n_u(x) – n_l(x); H_a(x) = H_u(x) + H_l(x); H_b(x) = H_u(x) – H_l(x). Обозначим пять уравнений системы (9) последовательно сверху вниз как E₁, E₂, E₃, E₄, E₅. Подставим в уравнения системы (9) выражение (10) для M(x) и перейдем к новой системе уравнений: E₁₀ = E₁ – E₂, E₂₀ = E₁ + E₂, E₃₀ = E₃ – E₄, E₄₀ = E₃ + E₄. Получим:

$\begin{array}{l}{E}_{10}:l^2{b}^{\text{'}\text{'}}-b+{H^{\text{'}}}_a=\mathrm{0,}\\ E_{20}:\left(l^2+\frac{A{h}^2}{4}\right)a^{\text{'}\text{'}}-a+\left(1+\frac{A}{2}\right){H^{\text{'}}}_b=\mathrm{0,}\\ E_{30}:\left(1+s\right){H^{\text{'}\text{'}}}_b-\left(\frac{A}{h^2}+\frac{k_s}{2}\right)H_b+\\ +\frac{k_s}{2}{H}_a-\left(1+\frac{A}{2}\right)a^{\text{'}}+\frac{2A}{h}=\mathrm{0,}\\ E_{40}:\left(1+s\right){H^{\text{'}\text{'}}}_a-\frac{k_s}{2}{H}_a+\frac{k_s}{2}{H}_b-b^{\text{'}}=0.\end{array}$ (11)

Из уравнения E₃₀ выразим H_a:

$\begin{array}{c}{H}_a=\frac{1}{k_s}\left(\left(2+A\right)a^{\text{'}\text{'}}-2\left(1+s\right){H^{\text{'}\text{'}}}_b+\right.\\ \left.+\left(\frac{2A}{h^2}+k_s\right){H^{\text{'}}}_b-\frac{4A}{h}\right),\end{array}$ (12)

и подставим в уравнения E₁₀, E₂₀, E₄₀. Получим уравнения E₁₁, E₂₁, E₄₁ соответственно:

$\begin{array}{l}{E}_{11}:2\left(1+s\right){H^{\text{'}\text{'}\text{'}}}_b-\left(\frac{2A}{h^2}+k_s\right){H^{\text{'}}}_b-\\ -\left(2+A\right)a^{\text{'}\text{'}}-l^2{k}_s{b}^{\text{'}\text{'}}+k_{s}b=\mathrm{0,}\\ E_{21}:\left(l^2+\frac{A{h}^2}{4}\right)a^{\text{'}\text{'}}-a+\left(1+\frac{A}{2}\right){H^{\text{'}}}_b=\mathrm{0,}\\ E_{41}:2{\left(1+s\right)}^2{H}_b^{\left(4\right)}-2\left(1+s\right)\left(\frac{A}{h^2}+k_s\right){H^{\text{'}\text{'}}}_b+\\ +\frac{A{k}_s}{h^2}{H}_b-\left(1+s\right)\left(2+A\right)a^{\text{'}\text{'}\text{'}}+\\ +\left(1+\frac{A}{2}\right)k_s{a}^{\text{'}}+k_s{b}^{\text{'}}-\frac{2A{k}_s}{h}=0.\end{array}$ (13)

Из уравнения E₂₁ выразим производную H_b′:

${H^{\text{'}}}_b=-\frac{4l^2+A{h}^2}{2\left(2+A\right)}{a}^{\text{'}\text{'}}+\frac{2}{2+A}a$, (14)

подставим H_b′ в уравнение E₁₁ и в продифференцированное по x уравнение E₄₁. Получим уравнения E₁₂ и E₄₂ соответственно:

$\begin{array}{l}{E}_{12}:-\frac{8A\left(l^2-h^2\right)+8h^{2}s+k_s{h}^2\left(4l^2+A{h}^2\right)}{2h^2{k}_s\left(2+A\right)}{a}^{\text{'}\text{'}}+\\ +\frac{\left(1+s\right)\left(4l^2+A{h}^2\right)}{k_s\left(2+A\right)}{a}^{\left(4\right)}+\\ +2\frac{2A+k_s{h}^2}{k_s{h}^2\left(2+A\right)}a+l^2{b}^{\text{'}\text{'}}-b=\mathrm{0,}\\ \\ E_{42}:-\frac{4h^{2}s+4A\left(l^2-h^2\right)+k_s{h}^2\left(4l^2+A{h}^2\right)}{h^2\left(1+s\right)}{a}^{\left(4\right)}+\\ +2\frac{2A\left(1+s\right)+k_s{h}^2\left(1+2s\right)+A{k}_s\left(l^2-h^2\right)}{h^2{\left(1+s\right)}^2}{a}^{\text{'}\text{'}}-\\ -\frac{2A{k}_s}{h^2{\left(1+s\right)}^2}a-\frac{k_s\left(2+A\right)}{{\left(1+s\right)}^2}{b}^{\text{'}\text{'}}+\left(4l^2+A{h}^2\right)a^{\left(6\right)}=0.\end{array}$ (15)

Из уравнения E₄₂ выразим вторую производную b″:

$\begin{array}{l}{b}^{\text{'}\text{'}}=\frac{{\left(1+s\right)}^2\left(4l^2+A{h}^2\right)}{k_s\left(2+A\right)}{a}^{\left(6\right)}-\frac{2A}{h^2\left(2+A\right)}a-\\ -\frac{\left(1+s\right)\left(4A\left(l^2-h^2\right)+k_s{h}^2\left(4l^2+A{h}^2\right)+4h^{2}s\right)}{h^2{k}_s\left(2+A\right)}{a}^{\left(4\right)}+\\ +2\frac{2A\left(1+s\right)+A{k}_s\left(l^2-h^2\right)+k_s{h}^2\left(1+2s\right)}{h^2{k}_s\left(2+A\right)}{a}^{\text{'}\text{'}},\end{array}$ (16)

и подставим в дважды продифференцированное по x уравнение E₁₂. В результате получится линейное однородное дифференциальное уравнение с постоянными коэффициентами восьмого порядка относительно функции a(x) вида

$Q_8{a}^{\left(8\right)}+Q_6{a}^{\left(6\right)}+Q_4{a}^{\left(4\right)}+Q_2{a}^{\text{'}\text{'}}+2Aa=\mathrm{0,}$ (17)

где

$\begin{array}{l}{Q}_8=\frac{l^2{h}^2{s}^2}{k_s}\left(4l^2+h^2\right)+\frac{l^2{h}^2}{k_s}\left(4l^2+A{h}^2\right)\left(1+2s\right),\\ Q_6=-4l^4{h}^2\left(s+1\right)-l^2{h}^4\left(s+A\right)-\frac{h^2{s}^2}{k_s}\left(8l^2+A{h}^2\right)-\\ -\frac{sA}{k_s}{\left(2l^2-h^2\right)}^2-\frac{8h^2{l}^{2}s}{k_s}+\frac{4A{l}^2}{k_s}\left(h^2-l^2\right),\\ Q_4=h^2\left(8l^2+A{h}^2\right)s+\frac{A}{2}{\left(2l^2-h^2\right)}^2+4h^2{l}^2+\\ +\frac{4h^2{s}^2}{k_s}+\frac{4sA}{k_s}\left(2l^2-h^2\right)+\frac{4A{l}^2}{k_s},\\ Q_2=2A\left(h^2-2l^2\right)-4s\left(h^2+\frac{A}{k_s}\right).\end{array}$ (18)

В уравнение (17) подставлялась функция a(x) вида a(x) = Ce^lx, где C и l – постоянные коэффициенты. Это приводило к характеристическому полиному уравнения (17) следующего вида:

$Q_8{\lambda }^8+Q_6{\lambda }^6+Q_4{\lambda }^4+Q_2{\lambda }^2+2A=0.$ (19)

Комплексные корни l₁, ..., l₈ характеристического полинома (19) находились аналитически, однако соответствующие выражения очень громоздки, поэтому здесь не приводятся. Поскольку в характеристический полином входят лишь четные степени l, то l₁ = –l₅, l₂ = –l₆, l₃ = –l₇, l₄ = –l₈; кроме того, пары корней l₁, l₂ и l₃, l₄ – комплексно сопряженные. Общее решение уравнения (17) записывалось в виде:

$\begin{array}{c}a\left(x\right)=C_1{e}^{{\lambda }_{1}x}+C_2{e}^{-{\lambda }_{1}x}+C_3{e}^{{\lambda }_{2}x}+C_4{e}^{-{\lambda }_{2}x}+\\ +C_5{e}^{{\lambda }_{3}x}+C_6{e}^{-{\lambda }_{3}x}+C_7{e}^{{\lambda }_{4}x}+C_8{e}^{-{\lambda }_{4}x},\end{array}$ (20)

где C₁, ..., C₈ – постоянные комплексные коэффициенты, которые следует определить из граничных условий. Далее a(x) подставлялось в уравнение (16), откуда после двух последовательных интегрирований по координате x получалось выражение для b(x); двойное интегрирование увеличивало число неопределенных коэффициентов до десяти. Выражение (20) также подставлялось в уравнение (14), откуда после интегрирования по координате x получалось выражение для H_b(x) и попутно возникал еще один (одиннадцатый) неопределенный коэффициент. Полученные выражения для H_b(x) и для a(x) далее подставлялись в уравнение (12), что позволило найти H_a(x). Затем определялись функции n_u(x), n_l(x), H_u(x), H_l(x) как n_u(x) = [a(x) + b(x)]/2, n_l(x) = [a(x) – b(x)]/2, H_u(x) = [H_a(x) + H_b(x)]/2, H_l(x) = [H_a(x) – H_b(x)]/2; они подставлялись в уравнение (10), из которого находилась функция M(x). Далее все найденные функции n_u(x), n_l(x), H_u(x), H_l(x), M(x) подставлялись в исходную систему уравнений Эйлера–Лагранжа (9). Из условий выполнения уравнений системы определялись дополнительные неопределенные коэффициенты, возникшие при интегрировании уравнений (14) и (16). Полученные таким образом в аналитическом виде функции n_u(x), n_l(x), H_u(x), H_l(x), M(x) являются общим решением системы уравнений Эйлера–Лагранжа (9).

Предполагалось, что деформации должны быть ограничены всюду и затухать при удалении от кромки поры (x → +∞). На неопределенные постоянные комплексные коэффициенты также налагались условия вещественности всех функций n_u(x), n_l(x), H_u(x), H_l(x), M(x) при любом вещественном аргументе x. Далее функции n_u(x), n_l(x), H_u(x), H_l(x), M(x) подставлялись в функционал энергии (8), и после интегрирования аналитически находилась энергия W_h.

Функционал упругой энергии мембраны (4) получен в предположении малости деформаций. В частности, предполагалось, что |n| << 1, |N| << 1, что означает, что малы отклонения директоров и нормалей от их направления в горизонтальном плоском невозмущенном бислое. Однако на экваторе поры липидные молекулы расположены практически горизонтально, и n = 1. Для формально корректного описания экваториальной области поры для нее возможно ввести другое референтное состояние, отличное от состояния плоского горизонтального бислоя. В качестве референтного состояния будем использовать плоский вертикальный монослой, в котором малы проекции директора и нормали на ось Oz, а не на ось Ox (рис. 1). Параметры, относящиеся к этому участку кромки поры, будем обозначать индексом «v»: проекции директора и нормали на ось Oz обозначим n_v, N_v соответственно; форму нейтральной поверхности вертикального монослоя будем задавать x-координатами ее точек, функцией H_v(z) (рис. 1). Для вертикального монослоя функционал упругой энергии записывается в виде:

$\begin{array}{c}{W}_v=\int dz\left\{\frac{k_c}{2}{\left(\frac{d{n}_v}{dz}+J_0\right)}^2-\frac{k_c}{2}{J}_0^2\right.+\\ +\left.\frac{k_t}{2}{\left(n_v+\frac{d{H}_v}{dz}\right)}^2+\frac{\sigma }{2}{\left(\frac{d{H}_v}{dz}\right)}^2\right\}.\end{array}$ (21)

Вариация этого функционала по функциям n_v(z) и H_v(z) приводит к двум уравнениям Эйлера–Лагранжа:

$\begin{array}{l}{l}^2\frac{d^2{n}_v}{d{z}^2}-n_v-\frac{d{H}_v}{dz}=\mathrm{0,}\\ \left(1+s\right)\frac{d^2{H}_v}{d{z}^2}+\frac{d{n}_v}{dz}=\mathrm{0,}\end{array}$ (22)

где l² = k_c/k_t; s = s/k_t. Проинтегрируем второе уравнение по z. Получим следующую систему:

$\begin{array}{l}{l}^2\frac{d^2{n}_v}{d{z}^2}-n_v-\frac{d{H}_v}{dz}=\mathrm{0,}\\ \left(1+s\right)\frac{d{H}_v}{dz}+n_v=d_0,\end{array}$ (23)

где d₀ – некоторый постоянный коэффициент. Из второго уравнения возможно выразить производную dH_v(z)/dz = (d₀ – n_v)/(1 + s) и подставить в первое уравнение системы (23). Получится изолированное линейное неоднородное уравнение с постоянными коэффициентами второго порядка на функцию n_v(z). Решение этого уравнения записывается в виде:

$\begin{array}{c}{n}_v\left(z\right)=d_2\exp \left[\sqrt{\frac{s}{1+s}}\frac{z}{l}\right]+\\ +d_1\exp \left[-\sqrt{\frac{s}{1+s}}\frac{z}{l}\right]-\frac{d_0}{s},\end{array}$ (24)

где d₁, d₂, – неизвестные действительные коэффициенты, которые следует определить из граничных условий. Подставив это решение во второе уравнение системы (23), найдем функцию H_v(z):

$\begin{array}{l}{H}_v\left(z\right)=\frac{l}{\sqrt{s\left(1+s\right)}}\left(d_1\exp \left[-\sqrt{\frac{s}{1+s}}\frac{z}{l}\right]-\right.\\ \left.-d_2\exp \left[\sqrt{\frac{s}{1+s}}\frac{z}{l}\right]\right)+\frac{k_s-s}{s\left(1+s\right)}{d}_{0}z+d_3,\end{array}$ (25)

где d₃ – неизвестный действительный коэффициент, который следует определить из граничных условий. Далее, подставляя полученные выражения для n_v(z), H_v(z) в исходную систему уравнений Эйлера–Лагранжа (22), получим, что первое уравнение системы (22) выполняется, только если d₀ = 0. Таким образом, общее решение уравнений Эйлера–Лагранжа (22) может быть записано в виде:

$\begin{array}{c}{n}_v\left(z\right)=d_2\exp \left[\sqrt{\frac{s}{1+s}}\frac{z}{l}\right]+d_1\exp \left[-\sqrt{\frac{s}{1+s}}\frac{z}{l}\right],\\ H_v\left(z\right)=\frac{l}{\sqrt{s\left(1+s\right)}}\left(d_1\exp \left[-\sqrt{\frac{s}{1+s}}\frac{z}{l}\right]-\right.\\ \left.-d_2\exp \left[\sqrt{\frac{s}{1+s}}\frac{z}{l}\right]\right)+d_3.\end{array}$ (26)

Решения (26) подставляются в функционал энергии (21), что после интегрирования позволяет найти упругую энергию вертикального монослойного участка.

Обозначим точку сопряжения вертикального монослойного участка с нижним монослоем горизонтального бислойного участка {X₁, Z₁}, а с верхним монослоем горизонтального бислойного участка – {X₂, Z₂} (рис. 1). Тогда интегрирование в (21) должно производиться по z в пределах от Z₁ до Z₂. В точках сопряжения налагались условия непрерывности нейтральных поверхностей и директоров вертикального монослойного и горизонтального бислойного участков. Остальные неопределенные коэффициенты находились из условия минимума полной энергии системы (W_h + W_v).

Можно произвольно зафиксировать координату X₁, поскольку это определяет положение мембраны как целого относительно начала выбранной декартовой системы координат. Остальные координаты сопряжения вертикального монослойного и горизонтального бислойного участка (X₂, Z₁, Z₂) находились путем численной минимизации полной энергии системы методом спуска по градиенту. Задавались некоторые стартовые значения координат ($X_2^0,Z_1^0,Z_2^0$). Затем на i-й итерации приближенно вычислялся вектор градиента упругой энергии системы в виде:

$\text{grad}\left(W\right)=\left[\begin{array}{l}\frac{W\left(X_2^{\left[i\right]}+{\delta }_0,Z_1^{\left[i\right]},Z_2^{\left[i\right]}\right)-W\left(X_2^{\left[i\right]},Z_1^{\left[i\right]},Z_2^{\left[i\right]}\right)}{{\delta }_0}\\ \frac{W\left(X_2^{\left[i\right]},Z_1^{\left[i\right]}+{\delta }_0,Z_2^{\left[i\right]}\right)-W\left(X_2^{\left[i\right]},Z_1^{\left[i\right]},Z_2^{\left[i\right]}\right)}{{\delta }_0}\\ \frac{W\left(X_2^{\left[i\right]},Z_1^{\left[i\right]},Z_2^{\left[i\right]}+{\delta }_0\right)-W\left(X_2^{\left[i\right]},Z_1^{\left[i\right]},Z_2^{\left[i\right]}\right)}{{\delta }_0}\end{array}\right],$ (27)

где d₀ – некоторое малое приращение (использовалась величина d₀ = 10^–10 нм); $X_2^{\left[i\right]},Z_1^{\left[i\right]},Z_2^{\left[i\right]}$ – значения координат X₂, Z₁, Z₂, соответственно, на i-й итерации. Значения координат на (i + 1)-й итерации находились согласно соотношению:

$\left[\begin{array}{l}{X}_2^{\left[i+1\right]}\\ Z_1^{\left[i+1\right]}\\ Z_2^{\left[i+1\right]}\end{array}\right]=\left[\begin{array}{l}{X}_2^{\left[i\right]}\\ Z_1^{\left[i\right]}\\ Z_2^{\left[i\right]}\end{array}\right]-{\delta }_C\text{grad}\left(W\right),$ (28)

где d_C – некоторая константа, выбираемая таким образом, чтобы обеспечить сходимость алгоритма; типичное значение d_C составляло ~2 × 10^–5 нм. Спуск по градиенту останавливался, когда каждая компонента вектора градиента по абсолютной величине оказывалась меньше 10^–5 k_BT/нм (k_BT ≈ 4 × 10^–21 Дж).

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Для получения количественных результатов использовались численные значения параметров, типичные для мембраны, сформированной из диолеоилфосфатидилхолина. Считалось, что модуль поперечного изгиба k_c = 10 k_BT [28]; спонтанная кривизна J₀ = –0.091 нм^–1 [34]; модуль наклона k_t = 40 мН/м [25, 29]; модуль латерального сжатия k_a = 133 мН/м [28]; толщина гидрофобной части монослоя h = 1.45 нм [24]; латеральное натяжение s = 0.1 мН/м. Величина K_s варьировалась в пределах от 0 до 10⁴ k_BT/нм⁴; в случае K_s = +∞ линейное натяжение было вычислено отдельно, с явно поставленным условием H_l(x) ≡ 0.

На рис. 2 показана форма мембраны вблизи кромки поры при значениях крутизны потенциала взаимодействия с подложкой K_s = 0 (рис. 2а), K_s = 10 k_BT/нм⁴ (рис. 2б), K_s = +∞ (рис. 2в). Величина K_s = 0 соответствует свободной мембране. В этом случае монослои мембраны находятся в одинаковых условиях, система симметрична, и межмонослойная поверхность плоская. При этом форма нейтральной поверхности нижнего монослоя существенно отклоняется от плоскости: расстояние от линии профиля поверхности до горизонтальной оси на рис. 2а немонотонно зависит от координаты x. При K_s = 10 k_BT/нм⁴ (рис. 2б) отклонение формы нейтральной поверхности нижнего монослоя от плоскости значительно меньшее; расстояние от линии профиля поверхности до горизонтальной оси на рис. 2б слабо зависит от координаты x. Форма межмонослойной поверхности отклоняется от плоскости, но достаточно слабо. Случай K_s = +∞ (рис. 2в) соответствует бесконечно сильной адгезии мембраны к подложке. Нейтральная поверхность нижнего монослоя становится плоской, а форма кромки поры существенно несимметричной.

 

Рис. 2. Рассчитанная форма нейтральных поверхностей монослоев и межмонослойной поверхности мембраны вблизи одномерной кромки поры. Кромка обладает трансляционной симметрией вдоль оси Oy, перпендикулярной плоскости рисунка xz. Положение начала координат выбрано таким образом, чтобы самой левой точке мембраны соответствовала координата x = 0. Формы рассчитаны при различных значениях крутизны потенциала взаимодействия с подложкой: а – K_s = 0 (свободная мембрана); б – K_s = 10 k_BT/нм⁴; в – K_s = +∞ (бесконечно сильная адгезия мембраны к подложке). На панели в стрелками показано рассчитанное направление директора n_u(x) в верхнем горизонтальном монослое, n_l(x) в нижнем горизонтальном монослое, n_v(z) в вертикальном монослое.

 

Зависимость линейного натяжения кромки поры от величины K_s показана на рис. 3а. В случае свободной мембраны (K_s = 0) рассчитанная величина линейного натяжения составляет g = 22.42 пН, что согласуется с экспериментально измеряемыми значениями [13]. При увеличении K_s линейное натяжение также монотонно возрастает (рис. 3), поскольку отклонение формы нейтральной поверхности нижнего монослоя от плоскости в этом случае требует затрат энергии. При K_s = +∞, что соответствует бесконечно сильной адгезии мембраны к подложке, рассчитанная величина линейного натяжения составляет g = 23.20 пН. Относительное изменение линейного натяжения при вариации K_s от 0 до +∞ оказывается меньше 3.5%. На рис. 3б показана зависимость линейного натяжения кромки поры от K_s в полулогарифмическом масштабе. Расчет проделан для величины K_s, варьируемой в пределах от 0 до 10⁴ k_BT/нм⁴. Линейное натяжение в случае K_s = +∞ вычислялось отдельно, с явно поставленным условием H_l(x) ≡ 0. Таким образом, разработанная теоретическая модель поры в мембране на твердой подложке предсказывает чрезвычайно слабое влияние взаимодействия с подложкой на величину линейного натяжения – основную энергетическую характеристику кромки поры. Результаты наших расчетов позволяют обосновать использование модельных мембран на подложках для исследования процесса образования пор различными методами, в первую очередь атомно-силовой микроскопии, и дают количественную оценку вызванного наличием подложки искажения измеряемого линейного натяжения кромки поры и производных от него величин.

 

Рис. 3. Зависимость линейного натяжения кромки поры от крутизны потенциала взаимодействия мембраны с подложкой: а – в линейном масштабе; б – в полулогарифмическом масштабе. K_s = 0 соответствует свободной мембране, K_s = +∞ – случаю бесконечно сильной адгезии мембраны к подложке. Относительное изменение линейного натяжения в полном диапазоне вариации K_s составляет менее 3.5%.

 

Конфликт интересов. Авторы декларируют отсутствие явных и потенциальных конфликтов интересов, связанных с публикацией настоящей статьи.

Источники финансирования. Работа выполнена при поддержке Российского научного фонда, грант № 22-15-00265 (П.В. Башкиров и М.В. Сумарокова) и государственного задания № 122030900065-6 (С.П. Дудик).

Соответствие принципам этики. Настоящая статья не содержит описания каких-либо исследований с участием людей или животных в качестве объектов.

About the authors

D. I. Kostina

National University of Science and Technology “MISiS”

Email: akimov@misis.ru
Russian Federation, Moscow, 119049

M. V. Sumarokova

Research Institute of System Biology and Medicine of Rospotrebnadzor

Email: akimov@misis.ru
Russian Federation, Moscow, 117246

S. P. Dudik

Research Institute of System Biology and Medicine of Rospotrebnadzor

Email: akimov@misis.ru
Russian Federation, Moscow, 117246

P. V. Bashkirov

Research Institute of System Biology and Medicine of Rospotrebnadzor

Email: akimov@misis.ru
Russian Federation, Moscow, 117246

S. A. Akimov

National University of Science and Technology “MISiS”

Author for correspondence.
Email: akimov@misis.ru
Russian Federation, Moscow, 119049

References

Supplementary files