附录：核心公式与理论推导

本文档是《皮肤屏障的”水之道”：分子模拟揭示脂质相共存如何稳定间质水》的技术附录，包含详细的公式推导、方法学细节和补充图表分析。

一、ABF（见上一篇）

二、渗透系数的计算方法详解

本文中，渗透系数（$k_p$）的计算基于非均匀溶解-扩散模型，并结合经验公式进行校准。

2.1 基于自由能和扩散系数的经典模型

理论上，渗透系数的倒数，即 resistance（$R$），可以通过对膜内各处的 local resistance 进行积分得到：

\[\frac{1}{k_p} = R = \int \frac{\exp(\Delta G(z) / k_B T)}{D(z)} \mathrm{d}z\]

这个公式的物理意义是，总的穿膜 resistance 是膜内每一点的 local resistance 之和。Local resistance 由两部分决定：

$\exp(\Delta G(z) / k_B T)$：这部分代表”溶解“的难度。$\Delta G(z)$ 是分子在膜内$z$位置相对于在水中的自由能（即PMF）。这个值越大，分子越不愿意待在这个位置，相当于溶解度越低，resistance 越大。
$1/D(z)$：这部分代表”扩散“的难度。$D(z)$ 是分子在膜内$z$位置的局部扩散系数。扩散越慢，resistance 越大。

2.2 本文采用的简化与经验校准模型

由于直接计算$D(z)$的复杂性和不确定性，作者采用了一种更巧妙的简化模型。他们发现，对于所研究的亲脂性小分子，其在膜内的平均扩散系数 $D$ 主要与分子量（MW）有关，且与经典的Potts-Guy经验公式（$D \sim \exp(-0.0061 \times \mathrm{MW})$）高度一致。因此，他们将渗透过程简化为由一个关键能垒控制的过程。

\[k_P = \frac{D}{\lambda_0} P_{liq}\]

2.3 公式的通俗解释

这个公式可以这样理解：一个分子的渗透系数 $k_P$ 由三个因素共同决定：

$D$（它能跑多快）：这是分子的平均扩散系数，主要由其大小决定。
$\lambda_0$（它要跑多远）：这是一个有效路径长度。它不只是膜的厚度，还考虑了分子在膜内迂回曲折的路径，因此通常比膜厚度大得多。这是一个需要通过实验数据来校准的经验参数。
$P_{liq}$（它进入”赛道”的概率）：这是最关键的创新点。作者假设，渗透并非在膜的任何地方都能发生，而是主要通过流动性更强的液态无序核心区。因此，$P_{liq}$ 代表了分子从有序区成功进入这个无序”赛道”的概率。这个概率可以通过分子穿过有序-无序界面所需的自由能垒 $\Delta G_{o/d}$ 来计算：

\[P_{liq} = \exp(-\Delta G_{o/d} / k_B T)\]

最终，作者通过对一系列已知渗透性的分子进行MD模拟，计算它们的 $\Delta G_{o/d}$ 和 $D$，然后与实验的 $k_P$ 值进行线性回归，最终拟合得到了经验参数 $\lambda_0 \approx 59 \mu m$，从而建立了一个完整的预测模型。

三、加热-退火模拟的详细过程

3.1 模拟的目的

加热-退火模拟是一种经典的计算方法，用于探索系统的亚稳态结构。在实验中，合成皮肤脂质样品时经常需要加热来加速脂质混合和相转变。因此，作者通过模拟这一过程来研究在高湿度条件下，水分子如何重新组织。

3.2 初始结构的构建

作者首先构建了一个由四层水合双层膜堆叠而成的大体系。具体步骤如下：

单个双层膜的准备：使用前面提到的1:1:2:2 CER[NS]/CER[EOS]/CHOL/FFA组成，构建一个平衡的水合双层膜（如图1所示的SPP模型）。
垂直堆叠：将这个双层膜沿着膜法线方向（Z轴）复制4次，形成4层双层膜的堆叠结构。每两层膜之间有一个水层分隔。
体系尺寸：
- 小体系：16 × 16 × 32 nm³
- 大体系：24 × 24 × 32 nm³
水/脂比：模拟了两种含水量：
- 5:1 水/脂比（较高湿度）
- 2:1 水/脂比（生理性湿度）
pH条件：模拟了两种pH：
- 低pH：所有游离脂肪酸（FFA）都质子化
- 中性pH：50%的FFA质子化，50%去质子化

3.3 加热阶段（95°C，0.25 μs）

温度升高：将体系从30°C升温至95°C。这个温度远高于大多数神经酰胺的熔点（通常在60-90°C）。
为什么选择95°C：
- 打破有序脂质链的堆积
- 增加脂质分子的动能和流动性
- 促进不同双层膜之间的接触和融合
- 加速水分子的重新分配
时间尺度：0.25 μs（250 ns）足够让脂质发生大规模重排，但不至于完全破坏膜结构。
观察到的现象：
- 相邻双层膜在多个接触点发生半融合（hemifusion）
- 外层脂质单层融合，但内层仍保持独立
- 原本分隔双层膜的水层被”挤压”

3.4 退火阶段（30°C，1.8 μs）

温度降低：将体系从95°C缓慢冷却回30°C（生理温度）。
为什么需要退火：
- 让系统从高温的无序状态”凝固”到某个亚稳态结构
- 观察水分子在冷却过程中如何重新组织
- 模拟实验中样品制备后的冷却过程
时间尺度：1.8 μs是一个相当长的弛豫时间，足以让脂质重新排列成稳定的构象。
最终结构：
- 高含水量（5:1）：形成连续的水通道，贯穿整个脂质基质
- 低含水量（2:1）：形成孤立的水滴，被脂质头基包裹在疏水核心中

3.5 为什么不形成标准的LPP

在退火后的结构中，虽然观察到了类似LPP的一些特征（如CER[EOS]的伸展构象），但整体上保留了显著的双层膜痕迹，没有完全转变为均一的13 nm厚的LPP结构。原因如下：

时间尺度限制：即使1.8 μs的模拟在计算上已经非常昂贵，但对于脂质的大规模重组（特别是长链神经酰胺的重排）来说，可能仍然太短。
缺乏层间模板：在真实的角质层中，角质细胞表面的共价结合脂质可能作为”模板”，引导脂质组装成LPP。模拟中缺少这种模板效应。
半融合是亚稳态：半融合状态本身就是一个能量局部极小值，系统可能”卡”在这个状态，需要更长时间或额外的驱动力才能进一步演化。

3.6 水滴与水通道的形成机制

关键洞察：在加热过程中，当相邻双层膜发生半融合时，原本位于膜间的水层被”困”在了融合的脂质核心中。退火后：

水含量高：水分子足够多，可以形成连续的柱状通道
水含量低：水分子被分散成多个孤立的球形水滴

这个结果表明，任何能引起膜局部结构剧烈重排的事件（热、机械应力、化学促渗剂）都可能将界面水”包裹”到疏水核心中，从而创造亲水性渗透路径。

四、水滴自由能模型的详细推导

4.1 模型的物理基础

文中提到：The free energy of the surface S of a water droplet was modeled as the sum of the interfacial tension with the lipid phase and the elastic bending energy of the surrounding lipid layer.

这个模型基于两个能量贡献：

界面张力能：水-脂质界面的存在需要能量（$\gamma$），类似于水滴在空气中的表面张力。
弯曲弹性能：包裹水滴的脂质头基需要弯曲，偏离其自然的曲率，这需要额外的能量。

4.2 完整的自由能公式

\[F(S) = \int_S \left[ \gamma + \frac{K_c}{2} (c - c_0)^2 \right] \mathrm{d}A_S\]

其中：

$\gamma$：水-脂界面张力（单位：mN/m 或 kcal/mol/nm²）
- 本文使用水-辛醇界面张力作为近似：$\gamma \approx 8.5 \pm 2$ mN/m
$K_c$：脂质的弯曲模量（单位：kcal/mol）
- 通过SPP双层膜的面积压缩模量计算：$K_A = 273 \pm 35$ mN/m
- 使用聚合物刷模型转换：$K_c = 9.5 \pm 1.2$ kcal/mol
$c = r_x^{-1} + r_y^{-1}$：总曲率，$r_x$ 和 $r_y$ 是两个主曲率半径
- 对于球形水滴：$c = 2/r$
- 对于圆柱形：$c = 1/r$（沿柱轴方向曲率为0）
$c_0 = r_0^{-1}$：脂质头基的自发曲率（spontaneous curvature）
- 从实验推导：$r_0 \approx 2.7$ nm
- 这是脂质头基在高湿度下”最舒服”的弯曲程度
$\mathrm{d}A_S$：表面积微元

4.3 球形水滴的自由能

对于半径为 $r$ 的球形水滴：

表面积：$S = 4\pi r^2$
曲率：$c = 2/r$

代入公式：

\[F(r) = 4\pi r^2 \left[ \gamma + \frac{K_c}{2} \left(\frac{2}{r} - \frac{1}{r_0}\right)^2 \right]\]

展开：

\[F(r) = 4\pi \gamma r^2 + 2\pi K_c \left[ 4 - \frac{4r}{r_0} + \frac{r^2}{r_0^2} \right]\]

4.4 寻找能量最小值

对 $r$ 求导并令其为零：

\[\frac{\mathrm{d}F}{\mathrm{d}r} = 8\pi \gamma r + 2\pi K_c \left[ -\frac{4}{r_0} + \frac{2r}{r_0^2} \right] = 0\]

整理得到最稳定半径 $r^*$ 满足：

\[\gamma r + \frac{K_c}{4} \left( \frac{r}{r_0^2} - \frac{2}{r_0} \right) = 0\]

代入数值（$\gamma = 8.5$ mN/m，$K_c = 9.5$ kcal/mol，$r_0 = 2.7$ nm），求解得：

\[r^* \approx 1.3 \text{ nm}\]

这与模拟中观察到的水滴平衡半径完美吻合！

4.5 物理意义

$r < r^*$：水滴太小，界面张力占主导，系统倾向于通过吸收更多水分子来增大半径，降低单位面积的界面能。
$r = r^*$：达到平衡，界面张力与弯曲能的竞争达到最优。
$r > r^*$：水滴过大，脂质头基被迫弯曲成比 $r_0$ 更大的曲率，弯曲能惩罚很大，系统倾向于释放水分子来缩小半径。

4.6 圆柱形通道的能量

对于半径 $r$、长度 $L$ 的圆柱形通道：

表面积：$S = 2\pi rL + 2\pi r^2$（侧面 + 两个端盖）
侧面曲率：$c = 1/r$
端盖曲率：$c = 2/r$

总自由能：

\[F(r, L) = 2\pi rL \left[ \gamma + \frac{K_c}{2}\left(\frac{1}{r} - \frac{1}{r_0}\right)^2 \right] + 2\pi r^2 \left[ \gamma + \frac{K_c}{2}\left(\frac{2}{r} - \frac{1}{r_0}\right)^2 \right]\]

4.7 形成通道的能垒

假设从一个 $r = r^* = 1.3$ nm 的球形水滴出发，保持半径不变，拉伸成长度 $L = 6$ nm（足以连接到邻近水滴）的圆柱：

初始能量：$F(r^*, L=0) \approx 0$（定义为参考点）
最终能量：$F(r^*, L=6 \text{ nm})$

计算得到：

\[\Delta F \approx 43 \text{ kcal/mol}\]

如果允许体积变化（即从周围吸收更多水），最优路径的能垒稍低：

\[\Delta F \approx 33 \text{ kcal/mol}\]

4.8 能垒的意义

稀有事件：在 $k_BT \approx 0.6$ kcal/mol（30°C）时，玻尔兹曼因子：

\[P \sim \exp(-33/0.6) \sim 10^{-24}\]

这意味着在平衡条件下，水通道形成是极其罕见的事件。

可促进性：但这个能垒不是不可逾越的。外部干预（如促渗剂、超声波、机械应力）可以提供额外的能量或降低能垒，显著提高通道形成的概率。

五、粗粒化力场参数化细节

5.1 SDK方法的9-6 Lennard-Jones参数

本研究中使用的粗粒化力场基于SDK（Shinoda-DeVane-Klein）方法，采用9-6 Lennard-Jones势能函数而非传统的12-6形式。这种选择能够更好地描述软物质体系的相互作用。

为了描述神经酰胺和游离脂肪酸的头基，作者从小分子的热力学数据（密度、表面张力、水合自由能）推导了新的力场参数：

核心参数表：

CG粒子类型1	CG粒子类型2	LJ ε（kcal/mol）	LJ σ（Å）
N（酰胺NH）	N	0.2430	4.0506
O（羰基C=O）	O	0.3233	3.7880
N	O	0.5393	3.6246
N	W（水）	0.9000	4.6100
O	W	0.6690	4.2166
COOH	COOH	0.6500	3.0000
COOH	W	0.7627	4.5418

其中：

N：酰胺NH基团（神经酰胺的鞘氨醇骨架）
O：羰基C=O（神经酰胺的酰胺键）
W：一个CG水粒子代表3个真实水分子
COOH：羧酸基团（游离脂肪酸头基）

5.2 参数化策略

小分子模型化合物：
- 甲酰胺（NH₂CHO）和N-甲基甲酰胺（CH₃NHCHO）：用于代表神经酰胺的酰胺头基
- 丁酸（CH₃(CH₂)₂COOH）：用于代表游离脂肪酸的羧基
拟合目标：
- 对角相互作用（同类型粒子）：拟合纯物质的密度和表面张力
- 与水的相互作用：拟合实验水合自由能
- 非对角相互作用：使用几何平均组合规则
pH条件限制：
- 由于CG水模型缺少偶极矩，无法稳定COO⁻等带电头基
- 所有CG模拟仅适用于低pH条件（FFA完全质子化）
开源资源：
- 完整力场参数：https://github.com/CG-it/ffdb-sdk
- 模拟输入文件生成工具（CG-it）：https://github.com/CG-it/CG-it
- 兼容LAMMPS软件

六、六方有序性分析

6.1 六方有序参数的定义

六方有序性（$

\psi_6

$）是描述脂质尾链在膜平面上二维排列规整程度的参数，定义为：

\[\psi_6 = \frac{1}{N_{neighbors}} \sum_{j=1}^{N_{neighbors}} e^{i6\theta_j}\]

其中 $\theta_j$ 是第 $j$ 个最近邻原子相对于中心原子的角度。

物理意义：

**$ \psi_6 = 1$**：完美的六方晶格（固态有序，gel相）
**$ \psi_6 = 0$**：完全无序（液态无序，liquid-disordered相）
**$0 < \psi_6 < 1$**：液-固共存或液晶相

6.2 胆固醇的流动化作用

通过计算不同胆固醇含量下的六方有序参数，揭示了胆固醇对脂质膜流动性的影响：

胆固醇含量	AA模拟	CG模拟	相态
0%	0.75	0.55	固态有序gel
30%	0.48	0.50	固-液共存
50%	0.42	0.40	液态无序Ld