蛋白模拟课笔记

一圈3.6个残基，螺距0.54nm i=4形成氢键

不易形成α螺旋：（1）酸性碱性集中电荷排斥（2）ILE异亮氨酸太大（3）Pro螺旋破坏者，因为氢键角度问题

二硫键需要检查一下，251kj/mol呢

＞=二次方的是短程力因为球形积分为常数

DH 理论推导泊松方程+玻尔兹曼分布 ->是个非线性方程，因此需要泰勒展开

Hydrophobic:

水在表面排成“规则“的水笼子，把一部分不自由的水分子挤出来，变成自由的水分子

结构->构象变化->蛋白-蛋白结合【【【【

debye-huckel 理论需要低浓度，需要低于0.1mol/L

24th Feb

自然界的奖励就是能量

miminal frustration 最小阻挫

energetic barrier 形成了不该形成的

;topological barrier ; 是native contact 形成了确实该形成的，但是形成顺序有问题,比如没装东西就合上了口袋

我好像想到了手性催组装的的讲法，很小的k1/k2就会产生显著效果

rugged floding funnnel: 要想成功折叠，必须漏斗足够深大于涨落，对native state 要有强的倾向

所以，一个疏水-亲水随意分布的（不考虑其他因素）序列就不会折叠

微观的参量能不能对应到宏观可测量量？

玻璃化温度，玻璃一样卡死了，s=0
$$
\frac{\delta E}{\Delta E}~ \frac{T_{fold}}{T_{玻璃化}}
$$

28th Feb

CASP 竞赛？怎么会有这么正义的比赛

学术界的主流方法：同源建模—+threading（都是有模板的）

最早的预测方法，swiss model

根据氨基酸序列对比数据库，根据序列相近程度算分

E-value: how a match is likely to arise by chance E> 1就是不可靠的，小于0.05就显著了

zhang server-I-TASSER

1）对齐Threading MUSTER/LOMETS 相当于前面的对比得到每一块的结构

评价用TMscore 为什么不用RMSD? 放在分母上，不会产生特别离谱的偏差

2）Assempbly 跑蒙特克罗模拟，把每个结构块串起来

3）然后做聚类，相似的结构多就是更合理的

对任意两个结构做比较，算RMSD？
在聚出来的2000个构象，把他们做平津，然后再跑一边MC

Contact 怎么算

两个方向的对角线分别为helix 和sheet

Co-evolution analysis 在进化过程中两个位点残基发生共同突变，说明它在三位结构里形成contact

“”

A-B, 如果A变弱了一点，B就需要变强一点

3rd Mar

Contact 预测

寻找直接相关direct coupling analysis DCA

即排除由于A-B, B-C导致的A-C

最大熵方法

$$
c\
S=-\Xi P(A_1, A_2 …)lnP(A_1, A_2 …)\
限定条件1：构造出来的概率必须和实际序列观察到的一致\
P_i(A_i)=f_i(A_i) —-h_i(A_i)\
P_{ij}(A_i,A_j)=f_{ij}(A_i,A_j) —-e_{ij}(A_i,A_j)\
限定条件2： \Sigma P(A_1, A_2 …)l=1
\构造拉格朗日乘子:
L=\
求偏导：\
如果定义整个系统的哈密顿量：
H=-\Sigma h_i(A_i)

\配分函数的一些性质：、、\
（1）\frac{\partial lnZ}{\partial h_i(A_i)}=\frac{1}{Z}\sum_{A_k,k\neq i} x_i\

如果\alpha=0. 没有关联，如果\alpha=1, 就是一开始的序列\
引入勒让德变换(把x,y 的函数变成m1,m2的）：
\lnZ–>G,
\=C_{ij}^{-1}
$$

dca.rice.edu/portal/dca/home

需要先构造MSA，在数据库中找到相似序列

蛋白的功能区几百万年基本不突变，比如actin, 这方法就没法用

如果两个位置一起突变，但是晶体结构却不一样，说明很可能体内是多聚体出现的（类似cc那种）

模拟函数中尽量偏向contact？

没有模板怎么办？切割成很短的20个氨基酸的 QUARK（和tasser一个入口）
alpha fold 第一步输出夹角分布（10度区间）+距离分布；距离分布可以对夹角分布做优化
- 在线版 colab.research.google.com/github/deepmind/alphafold

7th Mar

Normal mode analysis

解决的还是时间尺度问题，分解的还是不同振动模式

chem,rev,201,110,3,1463-1497

基本假设：在能量最低点‘

写出1. V(q) 2. 动力学方程
$$
F=ma=\frac{\partial V}{\partial q}
$$
V’’构建得到hessian 矩阵, 求解，频率就是矩阵的本征值，方向就是本征向量
$$
(M^{-1}H)U_k=w^2_k U_k
$$
一般选取一两个重要的模式，$\omega$越小越容易移动，所以越重要，对功能有用的，是全局的，低频的大幅度的震动

http://anm.csb.pitt.edu/

Enter the PDB id of your protein 输入1CLL

优点：并不需要高精度的晶体

缺点：1) 构象变化比较大的话，违背了弹性假设（cracking0) 2) 不知道一个关键运动如何线性叠加，3）不知道这些模式占比多少

Principle Component Analysis

拿轨迹文件构建hessian矩阵，再分析mode

eg, 设计药物靶点。这样如果蛋白两个地方振动有关，未必非要在结合RNA的地方设置靶点

10th Mar

蛋白和水分子要分开用力场

Verlet三个缺点速度和位置不同步位置和速度不耦合速度是二阶小量位置四阶，所以位置更精确（为什么要速度位置一致？控温算法时候改的是速度，不同步就麻烦了
蛙跳解决了精确度问题但是速度位置自然不同步，某一个时刻不可能同时知道速度和位置
velociy-verlet

pbc 盒子取多大？d>D/2 或r_cut<(D+2d)

r_cut=1.2nm d=1.0-1.5

蒙特卡洛模拟，可以看马达是什么事时候变到下一个态的，会到哪个态

14 Mar

KR DE并不是某个原子带了1个电荷，而是好几个原子都带partial charge
如果显性考虑了所有粒子才不需要debye
目前没有一个精确的立场能够描述DNA打开过程，计算机再好都没用用（然后老师面试博后就淘汰了）
1）所有第一代水分子力场都是刚性力场，区别仅仅在于电量的分配 SPC, TIP3P, TIP4P。 2）只要不是特别在意水，或被证实和水分子特别敏感，第一代就足够需要了 3）要和力场匹配，amber就用TIP3P

21 Mar

控温

berebdseb: 普通npt: 通过对速度进行纠正，并没有用正确的动能分布。
t推荐使用v-rescale
所以用了nosehooer,,严格的正则系综温度会有大的震荡
如果coupling太长，对热力学平衡态影响不大，但对动力学影响可能很大

在做控制温度时候水和蛋白分开考虑，蛋白质偏冷，水偏热（因为蛋白质自由度偏少）

24 Mar

2.控压？

通过调整什么？claussius virial theory-通过改变体积

Constraints in MD

时间间隔被你系统频率最高的运动限制，氢相连的化学键。但是氢又没什么用。。。。那就加个限制键长（C-H）等，然后只能提速一倍。。。2fs两个方法SHAKE LINCS

1)给初始力场和坐标2）确定力场

gmx energy -f runmd.edr 13-energy 16-tempereatuer

xmgrace energy,xvg

npt, nvt需要把蛋白固定,md只是没有这一行define

amber99sb ildn

含有核酸的话，需要bs0、bsc1修正

28 Mar

折叠到非折叠，势能升高了

热力学：求解自由能W=-kTlnP

平均力势：自由能在反应坐标上的投影
$$
W= -kTln\int e^{-\beta E}dq_1 dq_2 dq_3
\在q1 上的投影
W1= -kTln\int e^{-\beta E} dq_2 dq_3
\W1/W2=<\nabla E>
\-\nabla E=F(力是势能求导)
<\nabla E>=
$$
D.E.Shaw: 能用钱解决的就不是问题，只改进硬件。先去华尔街赚钱，然后暴力破解，每年少部分时间开放申请，你可以写个本子，选中了可以算7天2010 Science的最大意义：告诉你力场确实是准确的

伞采样

$$
E=E’+V, V=1/2k(Q-Q_0)^2， (k, Q_0自己选)
$$

“window”每个窗口加独立的二次函数，获得在那个窗口附近比较好的采样

为什么叫伞形？ -lgP
要求窗口之间要互相重叠
折叠到非折叠，势能升高了

热力学：求

要在contribut文件夹里

31 Mar

伞采样

每个地方的初始构象要接近k(d-d0)^2 的d0，这样加了简谐式就不会因为是能太大飞起

max_ndx -f  npt.gro
pull                     = umbrella
pull_ncoords            =1
pull_ngroups             = 2 
pull_group1_name              =r_1
pull_group1_name              =r_19
pull-geometry              = distance
pull-coord1-groups       =1 2 
pull_coord1_dim         = Y Y Y
pull_coord1_rate         = 0.01
pull_coord1_k                  = 1000

#伞采样的时候
pull_coord1_rate         = 0.00
pull_coord1_init         = 1.1;1.2,1.3...


-select 'com of r_1 and com of r_19'
gmx distance -f pull.xtc -s pull.tpr -n index.ndx

vi tpr.dat里面一系列的tor
vi f.datumb.xvg
wham -it  tpr.dat -if f.dat -min 0.9 -max 1.3

峰有重叠，就可以分析了，否则就失败了
可以调整 1）k,太小降低不了能垒，太大

http://www.mdtutorials.com/gmx/umbrella/05_pull.html

wham进行插值

2 Apr

增强采样：只能得到热力学信息，不包含热力学信息

1. Umbrella

2. Metadynamics

每隔一段时间加入高斯函数，把坑填平。问题是不知道什么时候就加够了，公式忘了记了

需要安装plumed->再安装gromacs->再打补丁

3.Replica exchange REMD副本交换

规则：细致平衡

X（X1,X2….)当前系统状态（坐标q，速度p）.交换坐标，调整速度（直接根据$\sqrt{T/T_0}p$

$\delta=[\beta_i -\beta_j][E(q_i)]-E(q_j)$ 构象i和构象j的势能差*系数
$$
p=1, \delta<0//

\exp(-\delta),\delta>0
$$
1)和伞形抽样相比，伞形是并行的，万一有一两个节点出了问题，就很烦。。。对计算资源要求很高(改进: single-copy tempering method)

2)由于大量水分子能量的涨落，掩盖了蛋白质本身构象的变化

也用wham进行插值

18 Apr

什么是好的反应坐标？

贝叶斯估计，看有没有很高的峰

究竟为什么粗粒度有效？

因为有水啊。。。因为“粗粒度”指的是力场的粗粒度

全原子力场来自于实验拟合+量化计算

Structure-based model…(Go model)

要保证native contact 比 nonnative contact 更稳定,才能抹平漏斗里的小坑。所以只给native加相互作用就行了

residue-based 又叫Cα model

我们建立这个模型的初衷不是研究化学性质，只是为了观察大幅度的构象变化，势能随便取

但只能消除energetic frusration 不能消除topological frustration

structure-based modelVS Target MD(直接把A拉到B)，：我并不关心末态，而是怎么跨垒的

Langevin ：老师是pieere 局里，居里夫人成为寡妇

涨落耗散定理

雨点打在车上这个例子真是太好了太好了

一些不知道什么时候的事情（好像是小组会)

jxy小组会

扩散系数作为变量，大就会两项完全分开，因为“有限时间”，超过0.66就会先发生别的转变
从chi可以导出表面张力系数sigma(arxiv)，分配系数也可推导chi
“动力学路径不同”
RPA 在平均场的基础上taylor展开到二阶，考虑到点-点关联，2体势
physical activity vs chemical activity
一一个脑洞：某个external field 周期性的？active crowding?
对流由表面张力引起
“这个图能回答提出的问题吗？”
反应扩散方程，如果出了自由能相加上external field “外在输入能力控制反应方向”
$$
\frac{\partial c}{\partial t}=f(c)+D\nabla^2c
f(c)+kc+kc^2
\frac{\partial F}{\partial c}=\mu(c)
$$
finite size scaling : 相互作用长度的3.5倍
03年PRL，行波，可激发excitable介质，+拓扑缺陷序参量：成功的次数
非平衡复杂体系，一定是各因素相互制约，只要坐标轴（序参量）选的好，一定是非单调的

（感觉生物的人肯定不知道）扩散和粘滞的退耦，平动和转动的decoupling
算每个东西都要反应一定的物理
Q: Cage 分布对。。。非均一因素。。。生物体系显然是非均已的，最好能找到提现非均一度的

21 Apr 又有实战啦

从PDB生成粗粒化力场（天哪我要是早两年遇到这个课就好了）

用网址，zizihttps://smog.rice.edu/cgi-bin/GenTopGro.pl
contact怎么算。 cutoff就是常规理解，shadow指的是两个离得近但是中间被别的挡住了
什么模型？all-atom 不是全原子，而是对力场进行了粗粒化

模拟

为什么不用能量最小，nvt npt?因为structer-based 本来就是能量最小的

integrator: sd 意思是启动朗之万模拟

他们竟然爱用LJ 10-12

-noddcheck: 告诉gromacs bond 长也没关系！

折叠温度：看Q，折叠态和非折叠态各占一半

VMD

protein 改成all

align

pdb->复制为rmp.pdb

top里保留pair部分并去掉第三列

1	`awk '{print $1,$2,$4,$5}' tem > BDD.cont`

不需要再看自由能曲线，直接得到歌态分布P，然后W=-KTlnP

抓住核心，native contact最重要，nonnative都不重要

25 Apr

决定马达蛋白ATP结合到前脚还是后脚：前脚开口要大一点，ATP进去了还会跑

vanilla: 原味冰激凌，所以是最早版本

Dual-basin model

为什么不用两个LJ相加？LJ的排斥会使得只能得到右边的basin 所以要用两个Gaussian 叠加

马达蛋白：力比较小的时候ATP释放速度慢，“贴封条”

testrun

top文件中的【pair】从原来的两个参数（epsilon sigma), 变成了一大堆高斯势参数(老师版的)，6表示只有一个basin有contact, 7表示

gmx analyze -f Q2.xvg -dist -bw 0.01

Sequence-dependent model

1, Miyazawa-Jernigan contact (M matrix )1）从蛋白数据库中，contact多的应该相互作用就强 2）quasi-chemical 近似，不考虑bond链接顺序，基本是散落的残基. 求最大熵(marcromolecules, 1985, 18, 534-552)

Skolnic contact (S matrix) 高斯链：头和尾的距离分布是高斯函数（polymer science, 1997,6,676-678)
Betancout-Thirumalai contac(B contact) 把thr -thr作为标准（0 ）

woc这不就是max entropy得到的疏水参数嘛。、。。

5 Mar

AWSEM

比residue-based 包含了多一点信息

Vtotal+Vcontact+Vburial+VHB+VFM+VDSB:

Vburial: 在内部还是外部参数不同

VHB help form helix sheet

VDSB:制造溶剂化空腔

VHB

能量涨落和漏斗深度对应宏观量folding temperature/glass tempetature >>1

contact太大，漏斗深但是δE也大；contact 太小，漏斗太浅

loop有时候会很重要和DNA结合（卡进去），类似disorder(然后作者让实验者mutate瞬间就结合不上去了)

实战

https://awsem.rice.edu/

右侧MSA：找同源序列过来用

右侧Coevolution: 之前的DC 也是数据库找’

http://frustratometer.qb.fcen.uba.ar/

configurational : 又考虑残基种类又考虑环境（距离） mutational只改变种类

绿线： minimal frustration 能量最低，稳定，进化成功

红线：功能区（功能要求不能稳定）或

下载加载：VMD, 打开pdb, extention->tk console -> file load -> tcl文件

Build your own CG model

关键是得到概率求自由能。V0=-kBTlng(r) 用自由能近似粗粒度模型的能量？

egTrpCage in urea

SB? cannot quantify urea effect 不包含化学因素

1）跑一组单个残基之间（比如一个Arg+ Lys) 计算g(r)

2）得到epslon~epsilon（urea)后用SB

JCP, 132, 17511(2010)

最大熵方法

U 前面必然是玻尔兹曼因子 λ1=β

按照正则系综的定义，

不断一定α 直到偏导为零，可以证明Γ是正定的，所以又唯一解

J..chem. Theory. Comput。 2012,8,10. 3445-3451

机器学习

acs central science 3019, 5, 755,767 jiang wang, clementi

必须保证什么东西不变（怎么写他的loss function)

要求全原子收到的力在粗粒化的投影，和他粗粒化模型势能求导的结果一样

9 May

蛋白-蛋白结合

结合过程中自由能先会上升，因为结合时候构象受到限制，导致熵的减小

有没有可能，IDP也有保守与步保守的区别，从数据库学出一些统计势场项

RosettaDock打分函数

coil和无规则卷曲常出现在活性位点

什么都不知道的盲目预测

clusPro.org

hdock.phys.hust.edu.cn

12 May

cluspro.org 先打分出比较好的1000个，再聚类“如果他是对的，一个个他是容易出现的”

如果认为疏水主导，hydrophobic favored

经典主导， elex…

如果什么都不知道，用balanced(有一些经验项)

不想用经验，就用vdw+elec

计算小分子药物结合，突变的影响：微扰

一系列含有不同比例 solvation free energy的哈密顿，把水分子之间的相互作用慢慢关掉

三种方法都是精确的，只是实际问题收敛速度不同

Q: 先关电场还是县官范德华？

电场，因为vdw还有体积排斥。。。先关就电厂发散了

chem.sci. 22016, 7, 207-218 aldeghi

建造一些热力学循环

May 16

尝试把A-B拉开，如果是拉的足够慢，应该有W=ΔA 否则因为有耗散，W>==ΔA

Jarzynski不等式

相空间， z=z(p,q)

孤立系统，烟花应该是决定性的，相空间演化轨迹不能相交
$$
\frac{\partial f }{\partial t}+\nabla j=0
\ \nabla j=\Sigma(\frac{\partial f}{\partial p} \dot p+\frac{\partial f}{\partial q} \dot q)
$$
刘维尔定理->f(z,t,s)不随时间变化

因为是孤立系统，W=首态和末态能量差
$$
\overline{exp(-\beta U)}=exp(-\beta \Delta A)
$$
但是指数的平均不是特别容易收敛

算两个蛋白的相互作用

不能有很多构象变换，否则diffusion control，算法失效

流体力学相互作用 1/r

May 19

随机过程

蛋白 3维跳跃和1维找寻的结合找到DNA的位点

Pn(M)(走n步后，位移在M的概率) ， 0.5时候概率分布最宽

联系生物过程，

往前跳概率为U, 往后跳概率为W，赌博发图第一个月
$$
p=ut, q=wt\
v=u-w=\frac{d}{dt}\

$$
需要含时的蒙特卡洛算法Gellasgie算法（不就是KMC吗）

第一个随机数用于确定下一步的时间，这种产生使得时间满足那个指数分布，第二个随机数选择事件。

Fisrt passage time

5.10 小组会 jxy

如果窄不容易错配

方程既包括离散项也包括连续项

波前运动，反应形成在边界上

化学波+自组装 active 增强自组装

一个脑洞：突变对noise的resist？

（科大学生喜欢把自己立于不败之地，基因里的挑刺）

熵是依赖于观察者的

2个droplet的距离作为反应坐标

融合过程受到熵势垒

（染色质环境（蛋白）对区室化的影响）一个提问题思路：把已有的问题倒过来

压力是单位体积的功

由于不交换，PV=nRT

大球只体现弹性不体现粘性

may 23

1. 平均首达时间

加入-1态后会对路径进行筛选，导致到达1的时间反而减少了

PMF vs 自由能

$$
平均力势是自由能在反应坐标上的投影\在q1 上的投影
W1, W1/W=<\nabla E>而力是势能求导(<\nabla E>=)所以叫平均力势
平均力势能 f(x,y)在x上的投影，平均的是y上的势能
$$

计算AB两态的自由能
全原子模型，用伞采样，在A和B 之间产生很多个初始构象，在不同的反应坐标区域内（有重叠）加上简谐势进行模拟，最后再用wham方法进行插值得到势垒
a) SB和直接用外力拉开有什么不同？
把A拉到B 是在关心初末态，拉伸的过程并非正常跨势垒，但structure-based 是想知道如何跨过势垒的

b)(宁这么问显然是不对)模型成功应当有其他实验结果对其进行验证，双态转换是input

transition 非常sharp的才叫two state
蛋白A大蛋白B小研究结合
选用微扰方法，小分子对大蛋白的哈密顿量可以看作是微扰;

假设尺寸可以做微扰，伞采样有中间态的概念，如果微扰只能得到自由能的差，如果系统比较小，围绕的精确度高于伞采样
研究IDP构象变化
metadynamics/replica exchange

5.12 大组会brc czy

szilard’s engine

每次测量都伴随熵

信息论，“我在黑板写板书，我写不需要力气，学生看也是不费力气，我擦掉就需要费力气”

“知识就是力量” knowledge is power (entropy is Work)

提问“我希望它什么样最好”

类似于高维的投影到低维距离变小

Maxwell demon 也是一种reset，只是确定性的

inducing: 原本不发生的区域诱导它能发生

前提，认为couple的机制不消耗。

都是想把无序涨落做成有用功

stochastic resonance： “机会留给有准备的人”

5.17小组会 wjx

transcription factor

”反应体系动力学不同，结论不同“

振荡增强化学反应。如果边界处速率最高，则震荡local就是一种相分离

（Q: 是否能查到需要边界处的反应？或者需要慢进行的反应？）

结论：多步反应、高级数反应：正反馈机制，非线性
（Q: 不在Rg-Tc线上的便宜原因是什么？）

直接关联函数是能量量纲->两体势+三体势。。。
R水合弱于K

软件基础，物理，前沿都有了，专题也很有趣

writing

蛋白模拟课笔记

http://home.ustc.edu.cn/~lingerli/2022/03/31/Notes Protein/

Author

Carolinge Li

Posted on

March 31, 2022

Licensed under

Literature records Previous

Hexo 使用指南 Next

蛋白模拟课笔记

24th Feb

28th Feb

最早的预测方法，swiss model

zhang server-I-TASSER

Contact 怎么算

3rd Mar

Contact 预测

最大熵方法

7th Mar

Normal mode analysis

Principle Component Analysis

10th Mar

14 Mar

21 Mar

24 Mar

28 Mar

伞采样

31 Mar

2 Apr

1. Umbrella

2. Metadynamics

3.Replica exchange REMD副本交换

18 Apr

什么是好的反应坐标？

究竟为什么粗粒度有效？

Structure-based model…(Go model)

涨落耗散定理

21 Apr 又有实战啦

从PDB生成粗粒化力场（天哪我要是早两年遇到这个课就好了）

模拟

VMD

25 Apr

Dual-basin model

testrun

Sequence-dependent model

5 Mar

AWSEM

实战

Build your own CG model

9 May

蛋白-蛋白结合

什么都不知道的盲目预测

12 May

计算小分子药物结合，突变的影响： 微扰

建造一些热力学循环

May 16

Jarzynski不等式

算两个蛋白的相互作用

May 19

随机过程

联系生物过程，

Fisrt passage time

5.10 小组会 jxy

（染色质环境（蛋白）对区室化的影响）一个提问题思路： 把已有的问题倒过来

may 23

1. 平均首达时间

PMF vs 自由能

计算AB两态的自由能

a) SB和直接用外力拉开有什么不同？

蛋白A大蛋白B小 研究结合

研究IDP构象变化

5.12 大组会brc czy

5.17小组会 wjx

计算小分子药物结合，突变的影响：微扰

（染色质环境（蛋白）对区室化的影响）一个提问题思路：把已有的问题倒过来

蛋白A大蛋白B小研究结合