CN112149328B - 一种用于模拟分子化学趋向运动的程序算法 - Google Patents

Abstract

本发明公开了一种用于模拟分子化学趋向运动的程序算法，通过建立有关受体与配体交叉扩散的有限元模型，从统计热力学基本框架出发，严谨推导出统摄分子趋化行为的解析方程；开发可推演该特殊运动进程的核心算法，以流场配置和流体性质为输入参数，对分子趋化系统给予数值模拟，用于预测并检验实验结果。本发明遵循“第一性原理”，从基本的统计热力学定理出发，完备地推导出描述分子趋化的支配方程，毋须试探性数学表达和经验性参数定义，建立起一般性模型；将运动方程转化为程序的核心算法，以实际体系的性质为基础，高效模拟出具体分子结合体系的运动趋势及有关效应，并输出定量预测结果，用于数据挖掘比对。

Description

一种用于模拟分子化学趋向运动的程序算法

技术领域

本发明属于复杂化学流体的有限元分析领域，尤其涉及一种用于模拟分子化学趋向运动的程序算法。

背景技术

微流控芯片中酶分子流向其底物的自发迁徙已经得以广泛报道，但该催化化学趋向运动的机制迄今悬而未决，导致该领域目前存在诸多争议性的观点，有的研究结论甚至相互对立。例如，2015年Nature报道了利用单分子荧光相关光谱，发现尿酶、漆酶、碱性磷酸酶在其本征的放热型催化过程中，自身扩散系数渐增，并与底物浓度、时间呈线性相关，拟合斜率与各酶对底物的转化频率(k_cat)有关，颇具指征性。但是，美国科学院院士SteveGranick在运用超分辨显微成像分别跟踪单个尿酶和漆酶时，观察到的却是离散酶个体在逐渐远离其对应底物，呈现反(逆)趋化现象，由此提出了酶的迁移率应受k_cat和米氏常数(K_M)共同制约。其后，哥伦比亚大学Henry Hess课题组通过动态光散射实验，补充道“酶在催化转化阶段的变构与亚基解离，也会不同程度地干预大分子的实时扩散；部分酶的扩散系数并不随底物浓度、催化时间而变化”。另一方面，虽然催化反应的焓变一般较为显著，但德国Max Planck研究所的Golestanian教授精确测量却发现，实际反应微区的温度仅提升了不足1/1000K；他还发现糖酵解代谢中的二磷酸果糖酶催化果糖-1,6-二磷酸这一吸热反应，亦具有放热型酶催化的全部趋化行为特征。-这两方面共同说明热效应不是主导分子趋化的主要因素。上述研究结论总体使酶趋化的成因分析愈发复杂。

为了阐明分子趋化的动力学成因，研发人员按照Michaelis-Menten模型和其准稳态近似，对酶催化历程中的两个基本阶段：(1)可逆的底物结合与(2)不可逆的产物转化，做了隔离、孤立考察，发现了离子和小分子通过非键且非特异性作用向较大物体接近的现象。比如，伊利诺伊大学香槟分校的Paul Braun团队通过在平面静电场内以牺牲离子的运动自由度为代价，观察到带负电荷的荧光素自发向局域化季铵盐构成的人造“焓黑洞”汇聚，并将该现象开发用于加速靶DNA的表面捕获。该模型最早可追溯至2005年布朗大学LambertFreund发现的病毒侵入、质粒转染等过程中细胞膜表面受体向胞吞位点的主动聚集现象。已应用于膜分离，如巧用疏水半透膜附近的低渗透压，吸引染料向其表面不可逆地牢固粘附。反之，当包围蛋白或外泌体囊泡表面Debye层内一、二价离子梯度足够大时，这些聚合物或微纳颗粒的机动性也会成百上千倍地提高，产生与胶体尺寸、离子种类相关的扩散泳输运，已应用于蛋白质液/液相分离、脂域区划、光免疫治疗(photoimmunotherapy)，乃至原油开采等。

以上系统全部是由界面静电场为主的广义化学结合力所统摄，都侧重于阐述限域环境中各性质对研究主体速率的综合促进作用。其数学表达各异，例如牛津大学的J.Agudo-Canalejo博士通过线性组合的形式，在流动性本构方程中添加泛定“泳”项；而华盛顿大学的J.Schurr等修正了费克扩散定律，从Kirkwood-Buff积分出发，在多物理场耦合分析软件COMSOL中自定义了一个“增强扩散”(见公式1)项：

其中，D_R是受体的自由扩散系数，D_XD为所谓的交叉扩散系数，c_R和c_L分别是受体、配体的浓度。如此，绕过基本的热/动力学前提假设和严谨推导，以求通过形式特殊的试探表达，经人为参数修正、迭代逼近，近似解析受体分子“流”(J_R)随

的分布。但是，诸如扩散泳、扩散渗等本质皆属非平衡游移，与化学平衡主导的酶-底物特异性识别不符；另一方面，也与酶-底物单位时间内按化学计量比的键合方式存在着显著的差异。说到底，正是由于对分子运动图景的本质缺乏认识，才试探性地引入诸如“增强扩散”的概念；以及形而上地类比宏观运动学的动量守恒，认为产物释出导致酶分子的弹道式反冲和构象的质心偏移，继而想象出现了溶剂的瞬间差分应力(所谓“化学声表面波”)等，如此这般假想了一系列微观力学模型，致使其预测结果皆与实际观测所得存在着较大的偏差。鉴于此，为了更加严谨而准确地预测分子化学趋向行为，既而用于解释、指导相关原理与实践应用，本专利提出了一个新的分子趋化运动数值模拟技术。

发明内容

本发明的目的在于针对现有技术的不足，提供一种用于模拟分子化学趋向运动的程序算法。

本发明的目的在于针对现有技术的不足：一种用于模拟分子化学趋向运动的程序算法，包括数学模型、理论公式、计算程序和输入输出配置；

所述的数学模型包括分子化学趋向运动的抽象场景和其中的初始与边界条件；

所述的理论公式包括基本假设、参数定义和支配方程组；

所述的计算程序包括数据结构、数值模拟算法及其核心；

所述的输入输出配置包括流场状况、运动对象的状态参数及数据统计方法。

进一步地，所述的分子化学趋向运动的抽象场景包括2～3个入口、1个出口的微通道反应器、集体迁徙和反向汇聚两种浓度梯度生成方式，和以各流体汇流时的流道横截面作为采样观测的时空起点。所述的初始与边界条件包括受体和配体的系列进样浓度，和这两种物质在初始时刻的空间分布，尤其是在两相界面与管壁边缘处的浓度值。

进一步地，所述的基本假设包括由受体、配体和受配体复合物构成的理想溶液，其是一个处于热力学平衡态的巨正则系综，以及用于描述大分子相互作用过程的McMillan-Mayer溶液理论。所述的参数定义指用于系统建模的统计热力学巨配分函数的完备自变量，包括化学势、浓度、空间位置、自由能多项式、结合常数和扩散系数等。所述的支配方程组为分别有关受体和配体的作用力与位置关系的Fokker-Planck方程，和其派生的流动性本构方程，这两种方程的解析表达式。

进一步地，所述的数据结构包括Fortran和Mathematica两种开发环境中的体系变量定义和其对应的数值类型，各变量之间的线性依存关系，时间标定与坐标区划方法。所述的数值模拟算法及其核心包括支配方程组的程序语言实现，基于链表遍历更新的值传递方式，循环判据，边界值补偿办法，和与实验结果的偏差校正措施。

进一步地，所述的流场状况包括微流体几何参数、线流速和流经时长。

进一步地，所述微流体几何参数包括长宽高。

进一步地，所述的运动对象的状态参数包括受体、配体和二者复合物在模拟前后与位置有关的浓度和扩散系数，以及受体-配体的缔合/解离常数。

进一步地，所述的扩散系数由单分子荧光相关光谱和核磁共振扩散序谱测量获得，所述的缔合/解离常数由紫外可见分光光度法或荧光光谱法测得。

进一步地，所述的数据统计方法是从采样终点的受体分布积分中评定趋化距离，再用受体的趋化迁移量与配体的浓度作图，拟合函数关系，提取键合信息。

进一步地，计算流程包括以下步骤：

步骤1：变量定义，新建输出文件，依次输入如下参数：三流道器件全宽，模拟时长，中央流道的受体浓度，侧流道的配体浓度，受体与配体的解离常数，受体、配体和受体/配体复合物的扩散系数；建立数组，标定时间步和时间增量，初始化浓度数据。

步骤2：常量因子定值；根据支配方程中扩散系数、平衡常数与物质浓度三者的解析关系，构建核心算法，建立循环；计算受体、配体的自扩散与交叉扩散效果，通道边界的物质量赋值，计算交叉扩散系数和单步反应速率常数；值替换与链表更新。

步骤3：向指定目录输出结果，包括输入的参数配置档案，受体、配体、受体/配体复合物的最终浓度分布，三者的浓度梯度，受体/配体复合物占受体的摩尔百分比，以及与位置相关的浓度、扩散系数的曲率，等等，用于与实验结果比对和二次数据挖掘。

本发明的有益效果是：

(1)遵循“第一性原理”，从基本的统计热力学定理出发，完备地推导出描述分子趋化的支配方程，毋须试探性数学表达和经验性参数定义，建立起一般性模型。

(2)将运动方程转化为程序的核心算法，以实际体系的性质为基础，高效模拟出具体分子结合体系的运动趋势及有关效应，并输出定量预测结果，用于数据挖掘比对。

附图说明

图1为受体在配体梯度中的化学趋向运动模型示意图。微通道包含3个入口和1个出口，受体(receptor)由中间通道引入，配体(ligand)由旁侧流道进入，构成诱导受体集群迁徙的配体梯度。实验用典型的流场参数如图所示，如体积流速是50μL/h。模拟起始于靠近入口的虚线框所指示的汇流位置，并终止于出口附近的虚线采样区域(region ofinterest，简称ROI)，生成受体、配体和其产物在顺流而下各位置截面(相当于统一位置不同时刻)的横向浓度分布。

图2为单分子荧光相关光谱测量受体扩散常数的实验示意图。(A)为含不同浓度(0，5、50、500nM，5、50、500μM)咪唑配体的10mM pH 7.0 PBS中500nM ZnTSPP受体荧光的S型衰减痕迹的归一化自相关函数拟合；(B)是计算所得对应咪唑配体含量(从左到右)：0，5、50、500nM和5、50、500μM的ZnTSPP受体扩散速率的直方图。

图3为核磁共振扩散序谱的测量实验示意图。(A)为5mM ZnTSPP受体的NMR-DOSY谱图；(B)为100mM咪唑配体的NMR-DOSY谱图；(C)为5mM ZnTSPP与100mM咪唑以不同比例分散于D₂O中的混合物的NMR-DOSY谱图；(D)为100mM ZnTSPP与5mM imidazole的NMR-DOSY谱图。顶部谱图：沿垂直线剖切的各组分¹H-NMR全谱。右坐标：由水平线指示的扩散系数。特征化学位移(δ)：(δ_ZnTSPP：8.18、8.30、8.96ppm；δ_imidazole：6.86、7.52ppm)。

图4为10μM ZnTSPP受体在从低到高成倍配制的咪唑配体中采用不同方法测得的谱图；(A)为采用分光光度法测得的谱图，其中曲线a和b分别是100M imidazole与空白缓冲(10mM pH 7.0PBS)的UV-Vis吸收光谱；(B)为(A)中λ_ab＝422nm处的峰强度拟合所得的标准曲线；(C)采用荧光滴定法测得的谱图；(D)为(C)中λ_em＝658nm处的峰强度拟合所得的标准曲线。

图5为10μM受体面向不同浓度配体的模拟归一化梯度分布图；从中间向左边依次为0、1、10、100nM，1、10、100μM，1、10、100mM，和1M，以在通道宽方向的位置为函数。

图6为受体在不同配体梯度中相对空白的峰移动与配体浓度的关系、其Langmuir等温拟合，以及与实验结果的对比图。

图7为模拟所得的趋向作用力评价示意图。

具体实施方式

下面结合具体实施例和附图对本发明做进一步详述。

本发明通过建立有关受体与配体交叉扩散的有限元模型，从统计热力学基本定律出发，严谨地推导出统摄分子趋化行为的解析方程；开发可推演该特殊运动进程的核心算法，以流场配置和流体性质为输入参数，对分子趋化系统给予数值模拟，用于预测并检验实验结果。

本发明用于模拟分子化学趋向运动的程序算法，包括数学模型、理论公式、计算程序，和输入输出配置；

本发明中，所述的数学模型包括分子化学趋向运动的抽象场景和其中的初始与边界条件。作为优选，所述的分子趋化运动抽象场景包括几何纹样如Ψ型的通道(三流道器件)，有3个入口、1个合流点、1个反应室和1个出口；两种受体趋化运动模式：集体迁徙(collective migration)和反向汇聚(inverted focusing)，取决于受体和配体的梯度场设置，分别由入口的成分安排决定(从左到右)：配体(ligand)/受体(receptor)/空白(buffer)(集体迁徙)，和受体/受体+配体/受体(反向汇聚)；所述配体/受体/空白表示配体、受体、空白缓释液从左到右依次流入三流道器件的三个入口；所述受体/受体+配体/受体表示受体、受体/配体混合物、受体从左到右依次流入三流道器件的三个入口。所述的初始与边界条件包括受体和配体的进样浓度，零时刻这两种物质在合流点截面的横向分布函数(Dirac型和高斯型)，和在器件的非理想物理边界处的浓度衰减(边缘效应)。所有流场与流体状态参数皆可观测并量化。

本发明中，所述的理论公式包括基本假设、参数定义和支配方程组。其中，基本假设基于McMillan-Mayer溶液理论，从统计热力学理论出发，描述大分子的相互作用过程；具体而言，将上述分子趋化运动建模进一步简化、假定成一个由受体、配体和受配体复合物(即产物)构成的理想溶液，并抽象概括为一个处于热力学平衡态的巨正则系综；所述的参数定义包括描述该系综的统计热力学巨配分函数本身和完备的自变量，包括化学势、浓度、结合常数等标量，位置和扩散系数矢量，和与体系本征性质相关的自由能(内能)多项式；所述的支配方程组为根据巨正则系综基本前提条件和巨配分函数状态参数关系所推导的，有关受体、配体的趋化驱动力的Fokker-Planck方程，其衍生的流动性本构方程以及物质迁移率的解析表达。

本发明中，所述的计算程序包括数据结构、数值模拟算法及其核心。所述的数据结构是指Fortran和Mathematica两种程序开发环境中与模型参数一一对应的变量定义和数值类型，包括受体、配体和产物浓度数组，三者的扩散系数、结合常数、通道宽度、流速和流时长度实数型变量，以及时间步长、数组编号等整数型变量；其他变量的设立主要用作计算媒介，涉及扩散系数、结合常数的转换和变量标定，用于提高算法执行效率。所述的数值模拟算法及其核心包括支配方程组的程序语言实现，基于链表遍历更新的值传递方式、循环判据、边界值补偿策略和与实验结果的偏差校正公式。

本发明中，所述的输入输出配置包括流场状况、运动对象的状态参数及数据统计方法。作为优选，所述的流场状况包括微流道宽度、流经时长(反应时间)；所述的运动对象的状态参数包括受体和配体的初始浓度和扩散系数，以及受体-配体的解离常数；所述的扩散系数由单分子荧光相关光谱(Fluorescence Correlation Spectroscopy，简称FCS)和核磁共振扩散序谱(¹H-Nuclear Magnetic Resonance Diffusion-Ordered Spectroscopy，简称NMR-DOSY)测量获得，所述的解离常数由紫外可见分光光度法测得；所述的数据统计方法有关受体向特定浓度配体的趋化位移值，依次包括浓度采集、背景扣除、组间平均、归一化、第一累积量力矩计算、与对照组做差等若干步骤。

本发明的具体实施方式中，所述的用于分子化学趋向运动模拟的程序算法，其运算流程包括以下步骤：

步骤1，变量定义，新建输出文件，依次输入如下参数：三流道器件全宽、模拟时长、中央流道的受体浓度、侧流道的配体浓度、受体与配体的解离常数以及受体、配体和受体/配体复合物的扩散系数；建立数组，标定时间步和时间增量，初始化浓度数据。

步骤2，常量因子定值；根据支配方程组中扩散系数、平衡常数与物质浓度三者的解析关系，构建核心算法、建立循环；计算受体、配体的自扩散与交叉扩散效果，通道边界的物质量赋值，计算交叉扩散系数和单步反应速率常数；值替换与链表更新。

步骤3，向指定目录输出结果，包括输入的参数配置档案，受体、配体、受体/配体复合物(产物)的最终浓度分布，三者的浓度梯度，产物占受体的摩尔百分比，以及与位置相关的浓度、扩散系数的曲率，用于与实验结果比对和二次数据挖掘。

如图1所示，我们推演了一个模型以测试是否纯缔合/解离能够引发受体-配体之间的定向运动。该模型是一个包含受体、配体流的微通道，其有3个入口和1个出口；通过改变通道内的流场参数，如物质浓度、几何尺寸和体积流速，来调节层流扩散梯度。模拟起始于靠近入口的虚线框所指示的汇流位置，并终止于出口附近的虚线采样区域。在下述各实施例中，受体皆由中间通道引入，配体则由旁侧流道进入，构成诱导受体集体迁徙的配体梯度场。

实施例1

本实施例提供一种在图1所示建模场景中，受体/配体缔合的流体动力学分析方法。具体步骤如下：

步骤1：遵循Hill编撰的《统计热力学导论》(An Introduction to StatisticalThermodynamics)和Schellman的《大分子结合》(Macromolecular Binding)论文中所详述的McMillan-Mayer溶液理论，我们分析一个空间区域的平衡态性质，它的体积为V，比分子尺度大、但在宏观上小，并保持恒温T。

除了面向环境开放而化学势μ_s不变的溶剂(为简化模型，目前仅设一种溶剂，但整个形式体系完全适用于多溶剂条件)以外，该体积对三种溶质物种的流体开放u＝{P₀,L₀,PL}，分别对应游离的受体、游离的配体，和受体-配体复合物。既而可得该巨正则系综的热力学势：

Ω(μ_s,μ_u)＝-k_BTlnΞ＝-(P_s+Π)V (2)

其中，Ξ表示经Legendre变换与Helmholtz势相关的巨配分函数；P_s＝P_s(μ_s)，代表V中给定化学势μ_s的纯溶剂的压强，也体现每个移动中的分子的平动自由度；Π＝Π(μ_s,μ_u)定义为在固定μ_s由于溶质的加入而产生的渗透压。

各溶质物种的密度，ρ_u，由热力学关系译为：

其中β＝1/k_BT；a_u是溶质活度，与溶质化学势有关，根据

ρ⁰作为参比密度被引入是为去量纲，也即ρ⁰＝1M；而

是在既定μ_s的理想溶剂中、溶质u在标准状态化学势：

其中，

是物种u的化学势，其在理想气体混合物中的密度ρ_u；vap＝ρ⁰；而

是溶剂化自由能，与将一个溶质分子从理想气体混合物中转移到一个同等密度ρ_u的理想溶液中相关。

步骤2：我们做两个假设。第一，溶质组成一个理想溶液，如此有

这意味着ρ_u＝a_u。

第二，溶质在与它们结合相互作用有关的平衡态：P₀+L₀＝PL，这意味着

至此表明：

其中，

是无量纲的结合常量。由此，3个a_u中只有两个是独立变量。所以，分别用支配全部受体和配体浓度的有效化学势μ_P和μ_L建模这三种溶质更加有用。另外，

而

这可在半巨正则系综中予以方便地证明：

然后根据

其中，结合自由能被表达成结合多项式ξ的形式：

其中，

和

特指ρ_P和ρ_L，而

同上。该多项式代表了一般化的内能表达，可用具体公式替代，比如Monod-Wyman-Changeux方程。

步骤3：Fokker-Planck支配方程。将巨正则系综与空间中位置R的各个小体积V关联，并假设在各如此小的体积单元内存在局域化的热力学平衡。然而，在更大尺度上，局部密度及相应的μ_P和μ_L或许会在空间中变化。

根据它们的化学势梯度，有效物种在R处将会受到一个净作用力：

这两个方程中的第一项是标准焓作用力，将扩散驱离高浓度区域；第二项来自结合自由能，驱赶分子驶向它们的自由结合伙伴。方程(14)和方程(15)强调分子化学趋向行为必然与化学反应相联系。

从这一点出发，假想系统处于一个过阻尼情境之中，并留意其平均速度和粘度相关的迁移率，对应的“流”可表达为：

第一项对应Fick定律，其中D_P和D_L分别是受体和配体整体浓度的有效扩散常数。原则上，这些扩散常数根据游离与结合分子比率的改变而随空间变化，但在此忽略该复杂情况。第二项，以坐标依存的非线性分式定义如下：

其形式上既不同与前文中的“增强扩散”项，后者源自扩散泳(diffusiophoresis)或扩散渗(diffusion-osmosis)，及其他复杂的非平衡泳流现象(如Belousov-Zhabotinsky振荡反应)；实质又是受体与配体为达成彼此缔合的化学平衡，而由ΔG_rxn所转化的热力学力。

从局部质量守恒推导出有关受体、配体空间密度演化的本构关系：

方程(16)至方程(21)描述了受体与配体缔合所导致的流体动力学定向迁移效应。

实施例2

本实施例提供一种程序算法，用于数值模拟图1所示微通道内分子动态耦合所致化学趋向行为。基于上述模型，在FORCE 2.0 FORTRAN Editor(免费下载链接http://force.lepsch.com)开发环境中编写和编译了代码，具体算法流程如下：

步骤1：平衡态近似。在位置R，底物与产物之间瞬时的缔合与解离是前提假设，意味着反应物的侧向扩散成为决速步，这保证了该数值模拟随时间流逝向着化学平衡态演进：

k_f·c_P(R)·c_L(R)＝k_b·c_PL(R) (22)

基于质量守恒，有c_P(R)＝c_P0(R)-c_PL(R)，其中c_P0(R)对应受体的初始浓度，k_f和k_b分别是正向和逆向反应速率。上式经取代和简化后，可得

其中，K_D＝k_b/k_f是受体-配体复合物的解离常数；反过来，缔合常数K_A＝1/K_D＝c_PL(R)/[c_P(R)·c_L(R)]。

步骤2：坐标依赖的扩散和迭代。按照遍历原理，三维稳态流可约简至一维情形下有关在确定位置x的受体在水平截面的浓度分布。方程(16)中的J_P(R)可被改写、简化：

此处，Φ＝c_PL(x)/c_P0(x)，是复合的受体分子的百分数。更重要的是

这表明D_Pb表观上既与复合受体分子的摩尔分数、又与受体vs.配体摩尔比相关。

为追踪演化轨迹，用J_P的空间导数(Fick第二定律)来执行迭代操作：

该方程阐明了遍及所有链表更新c_P的策略，它的刷新速度(即时间步长)通过扩散层厚度

标定。

配体的浓度梯度(dJ_L(x)/dx)也会受到同步微扰，配体流也遵守类似方程(26)的表达：

通过与方程(23)相减，引入循环增量(Δc_P)，用来实时刷新每个循环的受体(P)、配体(L)的消耗以及复合物(PL)的生成：

程序以一套完备的输入配置执行初始化，相关信息如表1所列。

表1：用于程序初始化的模拟所需参数配置，包括名称、数值和单位

输入参数名称	赋值及单位
		微通道全宽	0.0234cm
模拟总时长	3.0s
		受体分子在中央泳道的浓度	1.0×10<sup>-5</sup>M
配体分子在旁侧泳道的浓度	0、1.0×10<sup>-9～0</sup>M
		解离常数K<sub>D</sub>	2.18×10<sup>-6</sup>M
受体的扩散系数	3.95×10<sup>-6</sup>cm<sup>2</sup>·s<sup>-1</sup>
		配体的扩散系数	1.55×10<sup>-5</sup>cm<sup>2</sup>·s<sup>-1</sup>
复合物的扩散系数	3.98×10<sup>-6</sup>cm<sup>2</sup>·s<sup>-1</sup>

步骤3：对沿通道而下的受体和配体流采用一种“活塞流/平推流”(Plug-Flow)近似；同时，通过一维差分Smoluchowski方程，接入有关跨层流浓度分布演化的核心算法，即方程(24)至方程(28)，并翻译为如下Fortran语言，历经三个环节：

(1)变量定义，新建输出文件，依次输入如下参数：三流道全宽，模拟时长，中央流道的受体浓度，侧流道的配体浓度，受体与配体的解离常数，受体、配体和受体/配体复合物的扩散系数；建立数组，标定时间步和时间增量，初始化浓度数据。

(2)常量因子定值；根据支配方程中扩散系数、平衡常数与物质浓度三者的解析关系，构建核心算法，建立循环；计算受体、配体的自扩散与交叉扩散效果，通道边界的物质量赋值，计算交叉扩散系数和单步反应速率常数；值替换与链表更新。

(3)向指定目录输出结果，包括输入的参数配置档案，受体、配体、受体/配体复合物(产物)的最终浓度分布，三者的浓度梯度，产物占受体的摩尔百分比，以及与位置相关的浓度、扩散系数的曲率，等等，用于与实验结果比对和二次数据挖掘。

实施例3

本实施例提供一种用于模拟图1所示场景中配体缔合诱导受体趋化的优化算法，与之匹配的模拟程序在Wolfram Mathematica v10.2开发环境中生成，具体步骤如下：

步骤1：考虑一个含有三种稀释成分：游离受体(R₀)、游离配体(L₀)和受体-配体(RL)复合物的溶液。假设在溶液中的每个空间点，R，配体和受体都处于键合平衡，其平衡解离常数，K_D：

其中，

是一个无量纲的结合常数，而且c⁰是一个参考浓度。

定义

和

分别作为受体和配体的总浓度；定义μ_R和μ_L为对应的有效化学势，统摄c_R和c_L的平衡态分布。对各物质整体而言，α＝R或L，有效化学势或可表达为：

在该表达式中，第一项是游离物种在溶液中的标准状态的化学势，第二项来自平动熵的标准贡献，而末项则由于键合相互作用降低化学势。该成键贡献以一个结合多项式表示，其对受体和配体或可表达为：

两个有效的物种，α＝R和L，都将经历一个热力学力，驱动它们向化学势减少的区域运动：

对该力的第一项贡献驱使扩散远离高浓度区域，以最大化平动熵。第二项贡献驱赶分子驶向它们的自由结合对象，以最小化它们的结合自由能。

步骤2：假定溶质运动是在一个过阻尼情境之中，那么净流是

第一项对应费克定律(Fick’s Law)，其中D_R和D_L分别是全部受体和全体配体浓度的有效扩散常数。特别地，这些有效扩散系数依赖于游离与键合分子的比率。

第二项描述向着自由结合对象的流，含有一个反映平衡结合常数的表观“交叉扩散”常数：

当前推导的交叉扩散反映的其实是一个产生自平衡结合作用的热力学力，而不是那些更为复杂的非平衡态泳流或流体动力学现象。

质量的守恒意味着有效的浓度根据下式演化：

对每个有效物种，α＝R和L。

步骤3：常量赋值：流道全宽，模拟时长，中央流道的受体浓度，侧流道的配体浓度，受体与配体的解离常数，受体、配体和受体/配体复合物的扩散系数；变量定义，标定时间步和时间增量，常数因子定值；建立数组，初始化浓度档案。

步骤4：以支配方程(34)至方程(38)为计算内核，建立迭代循环；计算受体、配体的自扩散与交叉扩散效果，通道边界量赋值；计算扩散系数和反应速率常数，链表更新；输出结果，包括受体、配体、受体/配体复合物(产物)的最终浓度分布，三者的浓度梯度。

步骤5：程序调试。

(1)因为没有配体流入通道时的浓度分布，鉴于配体流入通道时的浓度分布难以观测，故假定各物种横跨对应通道的起始浓度是恒定的；为了数值稳定，略微平滑在通道入口边界处的浓度。

(2)鉴于该算法取决于三个参数：平衡解离常数K_D，和受体与配体的有效扩散系数D_R与D_L，对这些参数做了独立于溶液组成的近似；系统地变化这些参数，力求复现实验中代表受体分子时空分布的荧光曲线。

经过详尽测试，模拟在5～60ns的积分时间步长和999～9072个格点的网络中是稳健的。将不同配体浓度下，在下端采样ROI中模拟所得的游离与结合受体总浓度的侧向分布做图并汇总，如图5所示。

实施例4

本实施例提供一套用于算法执行和数值模拟的输入配置(如表1所列)，以及各参数的获取途径，具体步骤如下：

步骤1：为求在观测过程的反复平移、聚焦过程中精确定位并完全复位图1中的ROI，在1024×1024像素²的窗口(面元binning＝1)内，绘制一组参考线，表达为二维坐标系端点的形式：(x₁,y₁)→(x₂,y₂)，特征线段有：

[1]十字交叉轴：水平线(0,510)→(1022,510)，垂直线(510,0)→(510,1022)；

[2]通道外壁：垂直左壁(420,0)→(420,1022)，垂直右壁(600,0)→(600,1022)，左侧入口斜截(0,0)→(420,425)，右侧入口斜截(1022,0)→(600,425)；

[3]垂直三等分线：(480,0)→(480,1022)、(540,0)→(540,1022)；

[4]水平横截线：(420,0)→(600,0)、(420,100)→(600,100)、(420,200)→(600,200)、…、(420,900)→(600,900)、(420,1000)→(600,1000)；

[5]合流交点：(420,425)→(600,425)。

专一采样位置是：合流点，与(420,1000)→(600,1000)的水平截线。

在10x物镜中，单位像素长度等于1.3μm(光学和数字的放大因子皆为1)。微通道宽度为234μm(600-420＝180pixels)，其他几何尺寸包括：通道高度/深度/厚度100μm，横截面积：2.34×10^-2μm²。

步骤2：将近出口的引流距离排除，测得合流处到末端ROI的距离。流动注射泵向三个入口灌输的连续体积流速都是50μL/h(微量注射器的内径由游标卡尺测量，为1.3mm)，计算得平均线速度约等于1.8mm·s^-1，据此可得流经整个通道所需的时间。

更准确的是，记录目标荧光微球(2.0μm)顺流而下直线游移的连续录像(帧数：20ms)，校正流动时间。≥3s的模拟总时长(反应时间)足够产生可分辨的分布差异。

步骤3：通过单分子荧光相关光谱(Fluorescence Correlation Spectroscopy，简称FCS)测定受体(一种荧光活性锌卟啉，间-四(4-苯磺酸基)卟吩锌(II)，简称ZnTSPP))，在低浓度配体(咪唑，imidazole)中的扩散系数。

在以倒置荧光显微镜(IX-71,Olympus)为中心的平台上操作进行单分子扩散的测量。在该模块化的仪器中，样品在高分辨三维压电台(NanoView,Mad City Labs)上定位。从High-Q激光发生器(PicoTRAIN,40μW)产生一束532nm,80MHz,5.4ps的脉冲激光，导入Olympus60x/1.2N.A.水浸物镜，受体的荧光穿过二向分束器(Z520RDC-SP-POL,ChromaTechnology)并聚焦于50μm大小、0.22N.A.的光纤(Thorlabs)基面上，尔后导入HFAC-26前置光电倍增管，最终发送至SPC-630时间分辨单光子计数(Time-Correlated Single-Photon Counting，简称TCSPC)板(Becker&Hickl)。

如图2A所示，当受体进出衍射极限的高斯聚焦体积(Gaussian Focal Volume)时，它的荧光强度(I(t))按TCSPC的“先进先出”(First-In First-Out，简称FIFO)模式变化。随机扩散所致的读数涨落在下式中自相关：

并通过如下的多元三维扩散级数展开：

其中，i表示可别受体的序号，N_i代表它们在任一给定时间的观测体积内的平均数密度，τ指代穿越时间，w是结构因子(其值在4～8之间)。

根据Stokes-Einstein方程(41)，用含100nM Rhodamine B(已知D_{Rhodamine B}＝2.8×10^-6cm²·s^-1)的50μL去离子水校正w：

其中k_B·T是玻尔兹曼常数(Boltzmann’s constant)与绝对温度的乘积，η为粘度。用Levenberg-Marquardt最小二次回归分析拟合形式如方程(40)的单组分曲线，

是唯一可调的参数。

对于扩散常数D_i的受体分子i^th，它将花费一个特征扩散期

跨越一个短半轴r约500nm的聚焦椭圆。该关系可表达如下式，由图2A中的衰减部分提取

计算D_i：

根据方程(39)可知，τ_D与单分子的绝对亮度无关，只有ΔI(t)才对泛函G(τ)有贡献；而激发单重态的寿命比τ_D更短，不会干扰后者的准确度，所以拟合时不考虑单重态的贡献。

如图2B和表2所示，原位将5μL含500nM受体与50μL若干低浓度的配体混合，使受体在液滴中的数密度等同于聚焦体积内≤(10±2)个分子，该条件下D_ZnTSPP几乎不变。

表2：不同浓度咪唑配体的存在对ZnTSPP受体扩散性质的影响

步骤4：如图3所示，通过核磁共振扩散序谱(¹H-Nuclear Magnetic ResonanceDiffusion-Ordered Spectroscopy，简称NMR-DOSY)测定非荧光活性配体和高浓度配体中受体的扩散系数。

以重水(D₂O,99.9％,Cambridge Isotope Lab)为溶剂，配制等摩尔浓度的ZnTSPP和咪唑(1mM)溶液，不含咪唑但ZnTSPP浓度不变的重水溶液作为对照组。将500μL溶液移至5mm耐高压核磁管内，并拧盖密封。用配备了液氮冷却探头(5mm CPP BBO 500S1)的BrukerAV-III-HD谱分析仪(500.2MHz,298K)，轮流记录样品的¹H-NMR和DOSY谱图。为获得准确的信号积分，用含四阶多项式的函数校正基线；所有化学位移以ppm为单位并相对残留水峰(δ_w＝4.70)记录位置。用TopSpin v3.2处理数据。

在静止磁场中，一系列强度渐增的电磁脉冲给各分子内的指定核自旋的相位编码。经过一个扩散延时Δt后，在给定δ时期内施加一个解码梯度脉冲，后者不会翻转已扩散了的分子的相位。共振强度I的衰减速率正比于扩散率D和脉冲振幅G，如下所示：

在此，I₀是无梯度脉冲的信号强度，而γ是所测核的磁旋比(MagnetogyricRatio)。

NMR-DOSY测量所需参数基于一个径向涡流延时脉冲序列，伴随一对双极化梯度脉冲和两个扰相梯度。有关测量条件罗列如下：PULPROG＝ledbpgp2s(1d),Δt(D20)＝0.2s,δ(P30)＝1800μs,16slices,Time Domain Size(TD)＝16K,Sweep Width(SWH)＝3.46ppm,Acquisition Time(AQ)＝4.73s,Center of Spectrum(O1)＝7.83ppm,Number of DummyScan(DS)＝4,Number of Scan(NS)＝32,Relaxation Delay(D1)＝1s,Gradient RecoveryTime(D16)＝0.0002s,Spoil Gradient Pulse(P19)＝600s。沿z轴的脉冲梯度(GPZ6[％])从其最大强度(0.50T·m^-1)的2％到95％线性递增，期间对于每个Δt都按2％信号留数的标准去优化δ。为减少涡流，设置一个额外的Δt(D21)＝5ms。为避免任何来自场强上升过程的热效应所致温度涨落，以350L/h流通空气以保持设备内恒温；但不旋转NMR管，以防出现一阶甚至更高阶的边频带。

步骤5：如图4A和图4B所示，受体与配体的结合常数K_D，可以经由紫外可见分光光度法测定：

其中，θ＝(ε_A-ε_F)/(ε_B-ε_A)，ε_A是给定溶液的吸光率除以受体的含量，ε_B和ε_F分别是结合与游离受体的摩尔吸光系数。上式可转换成：

在Shimadzu UV-3600光谱系统中完成滴定操作，用其配备的软件UVProbe 2.6记录数据。通过绘制1/(ε_A-ε_F)与1/[c_L·(ε_B-ε_F)]之间的标准曲线，K_D值能从线性拟合的斜率中直接获得。

如图4C和图4D所示，受体与配体的结合常数K_D，亦可经由荧光滴定法获取。在Edinburgh FLS920光致发光谱仪中完成测定，用其配备的软件记录光谱强度数据。将不同配体浓度淬灭受体荧光强度的对应变化，用下述朗缪尔等温吸附函数(LangmuirIsotherm)拟合：

其中，ΔI是受体在不同配体浓度中相对空白对照组的荧光强度变化量，B对应最大淬灭幅度，K_D是配体-受体复合物生成的平衡解离常数，c_ligand是配体浓度。

实施例5

本实施例提供一种算法模拟结果的数学处理方法，具体步骤如下：

步骤1：确立衡量受体趋化幅度的指标，并予以计算。

将模拟条件表示为α，具体指代受体和配体的总浓度。根据图1模型所示，溶液沿z方向流经微通道。以通道内壁为起点，定义横向位置为x，0≤x≤L，L是通道宽度。总荧光强度

是x的函数，在以z为中心、宽度为dz的微元内测量。定义

为观测到的强度最小值：

定义归一化了的强度，I_α(x)：

I_α(x)作为归一化几率分布，有长度的倒数单位，并可直接与荧光素的数密度关联。考虑如下两个前提：

I.所测强度与荧光基团强度线性相关，根据

其中，b_α是背景信号，ρ_α(x)是位置x处的数密度，而k_α是各发光体对信号的贡献，k_α和b_α与x独立。

II.背景噪声可由所测强度的最小值估算：

或等价为

若前述两个前提条件满足，则

而且

其中

N_α＝∫₀ ^Ldxρ_α(x) (53)是所测荧光分子的总数，即处于z且宽度为dz的小面元。再者

因此，I_α(x)dx就是在dx间隙内有关x的发光体份数。

既然I_α(x)已被恰当归一，那么化学趋向位移，μ_α，就可被规范度量，其表达如下：

μ_α＝∫₀ ^LdxI_α(x)x (55)

该式称为第一力矩或第一累积数，特指荧光分子的平均位置。

通过比较两个不同实验的I_α(x)来定义“趋化移动量”，δμ。两次模拟中的受体浓度相同，但其中一个实验包含配体(α＝L，ligand)，而另一个是不含配体的对照组(α＝C，control)。这样，化学趋向移动可表示为：

δμ＝μ_L-μ_C＝∫₀ ^Ldx[I_L(x)-I_C(x)]x (56)如果受体的趋化方向朝着左侧泳道中的配体，那么δμ<0；故取其绝对值，|δμ|，使该“集体迁徙”算法不依赖于配体的布局。

步骤2：如图6所示，将不同配体浓度诱导受体发生趋化运动的移动数据点，用下述朗缪尔等温吸附函数(Langmuir Isotherm)拟合：

其中，δμ是受体在不同配体梯度中相对空白对照组的趋化位移量，B对应最大化学趋向平移(即幅度)，K_D是配体/受体复合物的平衡解离常数，c_ligand是配体浓度。

按照方程(57)拟合模拟信号，确定K_D＝65.2μM，该模拟结果落在表1和图4不同实验所测定的K_D＝34.5～74.4μM范围内。

步骤3：对方程(16)做适当的数学处理，将作用力(即浓度梯度)与流关联，将D_Pb假定为D_P的有限倍数。通过类比Fick第一定律，可得如下试探性表达：

其中第一项描述自发扩散，第二项贡献于化学趋向，速度ν_P＝F_PD_P/k_BT。c_L对F_P的作用可通过联立方程(16)和方程(58)阐明：

与方程(18)相比，括号内的第二项可被改写成

其中缔合常数K_A是K_D的倒数。通过不定积分，可得F_P对c_L的解析依赖，对每个预先给出的c_P，有

F_P是单分子化学趋向力的平均。如图5所示，按c_L/K_D的比例可被划分为三段：

(1)当c_L＜＜K_D时，F_P→0，对应不明显的化学趋向；

(2)当

时，F_P～k_BTln2，与图5中的拐点匹配；

(3)当c_L＞＞K_D时，受体被过量配体包围，逐渐达到它们的最大配位数。在此情形下，将k_B＝R/N_A代入方程(61)得：

上式说明，受体与配体的缔合过程所转化的能量来自ΔG_rxn，有效的化学趋向要求ΔG_rxn＞RTlnc_L，即需要ΔG_rxn的转化以抗衡随机扩散的熵变。

通过方程(62)，可进一步计算出模拟系统的化学趋向力(功)，如图7所示。

Claims (8)

1.一种用于模拟分子化学趋向运动的方法，其特征在于，包括数学模型、理论公式、计算程序和输入输出配置；

所述的数学模型包括分子化学趋向运动的抽象场景和其中的初始与边界条件；

所述的理论公式包括基本假设、参数定义和支配方程组；

所述的计算程序包括数据结构、数值模拟算法及其核心；

所述的输入输出配置包括流场状况、运动对象的状态参数及数据统计方法；

所述的分子化学趋向运动的抽象场景包括2~3个入口、1个出口的微通道反应器、集体迁徙和反向汇聚两种浓度梯度生成方式，和以各流体汇流时的流道横截面作为采样观测的时空起点；所述的初始与边界条件包括受体和配体的系列进样浓度，和这两种物质在初始时刻的空间分布；

所述的基本假设包括由受体、配体和受配体复合物构成的理想溶液，其是一个处于热力学平衡态的巨正则系综，以及用于描述大分子相互作用过程的McMillan−Mayer溶液理论；所述的参数定义指用于系统建模的统计热力学巨配分函数的完备自变量，包括化学势、浓度、空间位置、自由能多项式、结合常数和扩散系数；所述的支配方程组为分别有关受体和配体的作用力与位置关系的Fokker−Planck方程，和其派生的流动性本构方程，这两种方程的解析表达式。

2.根据权利要求1所述用于模拟分子化学趋向运动的方法，其特征在于，所述的数据结构包括Fortran和Mathematica两种开发环境中的体系变量定义和其对应的数值类型，各变量之间的线性依存关系，时间标定与坐标区划方法；所述的数值模拟算法及其核心包括支配方程组的程序语言实现，基于链表遍历更新的值传递方式，循环判据，边界值补偿办法，和与实验结果的偏差校正措施。

3.根据权利要求1所述用于模拟分子化学趋向运动的方法，其特征在于，所述的流场状况包括微流体几何参数、线流速和流经时长。

4.根据权利要求3所述用于模拟分子化学趋向运动的方法，其特征在于，所述微流体几何参数包括长宽高。

5.根据权利要求1所述用于模拟分子化学趋向运动的方法，其特征在于，所述的运动对象的状态参数包括受体、配体和二者复合物在模拟前后与位置有关的浓度和扩散系数，以及受体−配体的缔合/解离常数。

6.根据权利要求5所述用于模拟分子化学趋向运动的方法，其特征在于，所述的扩散系数由单分子荧光相关光谱和核磁共振扩散序谱测量获得，所述的缔合/解离常数由紫外可见分光光度法或荧光光谱法测得。

7.根据权利要求1所述用于模拟分子化学趋向运动的方法，其特征在于，所述的数据统计方法是从采样终点的受体分布积分中评定趋化距离，再用受体的趋化迁移量与配体的浓度作图，拟合函数关系，提取键合信息。

8.根据权利要求1至7任一所述用于模拟分子化学趋向运动的方法，其特征在于，计算流程包括以下步骤：

步骤1：变量定义，新建输出文件，依次输入如下参数：三流道器件全宽，模拟时长，中央流道的受体浓度，侧流道的配体浓度，受体与配体的解离常数，受体、配体和受体/配体复合物的扩散系数；建立数组，标定时间步和时间增量，初始化浓度数据；

步骤2：常量因子定值；根据支配方程中扩散系数、平衡常数与物质浓度三者的解析关系，构建核心算法，建立循环；计算受体、配体的自扩散与交叉扩散效果，通道边界的物质量赋值，计算交叉扩散系数和单步反应速率常数；值替换与链表更新；

步骤3：向指定目录输出结果，包括输入的参数配置档案，受体、配体、受体/配体复合物的最终浓度分布，三者的浓度梯度，受体/配体复合物占受体的摩尔百分比，以及与位置相关的浓度、扩散系数的曲率，用于与实验结果比对和二次数据挖掘。

Images