斐波那契晶格生成

通过斐波那契晶格生成算法，我们可以在一个任意半径的球体表面进行均匀采样，得到给定的N个坐标点。基于这个晶格生成算法，可以对整个分子的所有原子做斐波那契晶格生成，再剔除共用的格点，即可得到分子表面（SASA）的均匀格点分布。

技术背景

在分子动力学模拟领域，蛋白质表面的一些几何特性，是非常重要的指标。典型的例如蛋白表面的结合口袋等，往往都会有内陷的几何特征。而如何构建出一个蛋白的溶剂可及表面（Solvent Accessible Surface Area，溶剂可及表面面积），是一个重要的研究方向。

有一个思路是，对蛋白体系中的每一个原子，构建一个半径为\(R_{van}+R_{sol}\)的球体，其中\(R_{van}\)为范德华半径，\(R_{sol}\)表示溶剂探针半径。那么这个构建出来的球体的最外层，如果不包含在其他原子的球体内部，就是对应的所谓溶剂可及表面（SASA）。这就是探针法计算SASA的一个思路，那么其中比较关键的一环，是给定范德华半径和溶剂探针的半径之后，如何在一个球体的表面做均匀采样。

大致方案有两种，一个是直接在球体表面做随机撒点，那么只要撒的点足够多，最后就可以得到一个密度均匀的球体表面的格点。还有一个当下最常用的算法——构建一个斐波那契晶格。这个算法，利用了无理数的特性，在球体表面构建出一条无限不循环的路径，然后根据特定的高度和角度，在这个路径上进行采样，就天然的得到均匀分布的球体表面样本点。而且还可以根据所需点数做非常灵活的调整。

斐波那契晶格

假如是在一个二维平面的圆周上，任意选择一个起点，将角度步长设置为一个\(2\pi\)的无理数倍，这样在圆周上行进的次数足够多的情况下，由于无理数的稠密性，就会在圆周上形成距离接近于均匀分布的点阵。常规的来说，用任意无理数都可以。但是不同的无理数，在收敛速度上会略有差别。根据一些论证，黄金分割数\(\Phi=\frac{1+\sqrt{5}}{2}\)会是一个比较好的选择，在比较少的点数下，也可以尽快的达到接近于均匀分布的密度。对应的角度为：

\[\phi=2\pi(1-\frac{1}{\Phi})=\pi(3-\sqrt{5})\approx 2.39996 rad\approx 137.5^\circ \]

那么得到，第\(k\)个点所在的角度为：

\[\theta_k=k\phi\mod 2\pi, k=0,1,2,3,...,N-1 \]

这些点永远不会周期性的重合。可以看一个少样本数的例子：

理论上来说，有50个点，就可以达到圆周上比较均匀分布的效果，但实际上跑起来：

emmm...其实效果也没有想象中的那么好，但是目前比较主流的算法就是这个斐波那契晶格生成算法了。

三维空间代码实现

前面看了两个二维的例子，实际上落实到三维空间的时候，还需要引入\(z\)这个维度：

\[\begin{aligned} z_k &= 1-\frac{2k+1}{N},k=0,1,2,3,...,N-1 \\ r_k &= \sqrt{1-z_k^2} \\ \theta_k &= k\phi\\ x_k &= r_k\cos\theta_k \\ y_k &= r_k\sin\theta_k \\ \end{aligned} \]

这个写法的物理图像，就是\(z\)这一轴从\(1\)取到\(\frac{1}{N}-1\)的均匀高度，同时使用前面提到的黄金角在\(x-y\)方向上进行旋转。对应的代码实现为：

import numpy as np

def save_to_xyz(path, xyz_name, atom_name='C'):
    atoms = path.shape[0]
    atom_name = [atom_name] * atoms
    with open(xyz_name, 'w') as xyz:
        xyz.write('{}\n'.format(atoms))
        xyz.write('{}\n'.format(xyz_name.replace('.xyz', '')))
        for i in range(atoms):
            xyz.write('{}\t{}\t{}\t{}\n'.format(atom_name[i], path[i][0], path[i][1], path[i][2]))

def fibonacci_sphere(n):
    golden_angle = np.pi * (3 - np.sqrt(5))
    i = np.arange(n)
    z = 1 - (2 * i + 1) / n
    r = np.sqrt(1 - z * z)
    theta = golden_angle * i
    x = r * np.cos(theta)
    y = r * np.sin(theta)
    return np.stack([x, y, z], axis=1)

if __name__ == "__main__":
    n = 4
    radius = 3
    xyz_name = "fib.xyz"
    nodes = fibonacci_sphere(n) * radius
    save_to_xyz(nodes, xyz_name)

这里我是将输出保存为xyz分子构象文件的格式了，然后用molstar这样的在线平台进行可视化：

在4个点的情况下，这个形状就很接近于正四面体。那么如果增加点数，整体看起来就会像是在球体表面均匀分布的一系列点阵：

最后把这个点阵信息输出给下一步的筛选模块进行处理，就可以得到最终的SASA表面。

总结概要

通过斐波那契晶格生成算法，我们可以在一个任意半径的球体表面进行均匀采样，得到给定的N个坐标点。基于这个晶格生成算法，可以对整个分子的所有原子做斐波那契晶格生成，再剔除共用的格点，即可得到分子表面（SASA）的均匀格点分布。

版权声明

本文首发链接为：https://dechinphy.github.io/posts/b185330b.html

作者ID：DechinPhy

更多原著文章请参考：https://dechinphy.github.io/