在完成docking之后，看著interaction的結(jié)果總感覺咱們的配體不可能有這么好的對接結(jié)果，所以分子動力學(xué)的驗證就應(yīng)允而生了。不僅能讓自己拿這個物質(zhì)準(zhǔn)備做實驗的心安下來，也能讓自己在驗證完這個藥物有效之后，進一步思考它的機理。今天就給大家分享一種open source的分子動力學(xué)方法，示例代碼加自己的蛋白+分子的實戰(zhàn)。

今天的教程來自teachopencadd，是一群大佬手把手教你怎么做docking，如果大家感興趣可以去https://projects.volkamerlab.org/teachopencadd/talktorials/T019_md_simulation.html閱讀原文；這篇內(nèi)容在教程的基礎(chǔ)上添加了如何準(zhǔn)備自己的配體和蛋白。

分子動力學(xué)模擬由三個分子組成的超分子在受限納米孔中的旋轉(zhuǎn)

分子動力學(xué)原理

如果要總結(jié)分子動力學(xué)在docking中的應(yīng)用，說人話就是：通過模擬微觀水平的分子之間的作用力，推斷出我們的配體在其中所有可能的存在情況。進一步判斷結(jié)合的穩(wěn)定性。 經(jīng)典分子動力學(xué) (MD) 模擬可用于研究涉及蛋白質(zhì)、脂質(zhì)、核酸、碳水化合物和小有機分子的過程。

汽車碰撞實驗是不是有點像分子動力學(xué)的物理版

那需要如何精確地描述分子的作用力呢？就涉及到其中幾個重要概念了：

力場

描述分子內(nèi)部和分子之間的原子間作用力。它們是參數(shù)方程，包含不同力（鍵伸縮力、范德華力等）的分量。參數(shù)值通常通過實驗推導(dǎo)，并根據(jù)所涉及的分子和模擬設(shè)置，在每個MD場景下發(fā)生變化。其結(jié)果是系統(tǒng)能量圖的數(shù)學(xué)描述，其中作用于每個粒子的力源于勢能相對于原子坐標(biāo)的梯度。

四大類 (生物) 分子力場 AMBER、CHARMM、GROMOS 和 OPLS ，可以參考這篇文獻：https://pubs.acs.org/doi/10.1021/acs.jcim.8b00042，詳細(xì)地描述了用于蛋白質(zhì)、核酸、脂質(zhì)、碳水化合物和小有機分子的 MD 模擬的固定電荷原子力場列表：

在這里我們將使用 AMBER 力場家族中的一個成員，該家族廣泛用于蛋白質(zhì)的 MD 模擬。它們的函數(shù)形式為：

是的，和大家一樣，這里我的大腦已經(jīng)停止運轉(zhuǎn)了。但幸好我們有可視化的圖片，從這里我們可以理解，通過在適應(yīng)的范圍內(nèi)改變分子的化學(xué)鍵的長度、角度和分子鍵的力度，同時考慮到分子間的力，就有機會模擬整個分子在環(huán)境中出現(xiàn)的狀態(tài)。

邊界條件

通常，分子系統(tǒng)是在充滿溶劑（例如水）的盒子中模擬的。這些盒子的尺寸有限，這會導(dǎo)致盒子邊界處或附近的分子出現(xiàn)問題。這些分子應(yīng)該與哪些分子相互作用？周期性邊界條件可以通過模擬理論上無限的系統(tǒng)來避免這種邊界偽影。因此，盒子一側(cè)邊界處的分子會與盒子另一側(cè)邊界處的分子相互作用。這模擬了模擬盒子被其自身復(fù)制品包圍的情況。在可視化此類MD模擬時，通?？梢杂^察到粒子從一側(cè)離開盒子（圖3）。然而，它們會同時以相同的速度重新出現(xiàn)在盒子的另一側(cè)。對于周期性邊界條件下的模擬，建議使用足夠大的模擬盒子，以便模擬的大分子不會與自身相鄰的圖像發(fā)生接觸。

具有周期性邊界條件的水分子的分子動力學(xué)模擬

MD模擬和藥物設(shè)計

MD 模擬為了解配體與其靶標(biāo)結(jié)合的高度動態(tài)過程提供了寶貴的見解。當(dāng)配體（或藥物）接近溶液中的大分子（蛋白質(zhì)）時，它會遇到一個不斷運動的結(jié)構(gòu)。此外，配體可能會誘導(dǎo)大分子發(fā)生構(gòu)象變化，以便最好地適應(yīng)小分子。靜態(tài)方法可能無法發(fā)現(xiàn)此類構(gòu)象。因此，在靜態(tài)無配體結(jié)構(gòu)中觀察到的結(jié)合位點，但可以通過MD模擬發(fā)現(xiàn)，有時被稱為隱蔽結(jié)合位點。使用 MD 模擬識別此類結(jié)合位點可以啟動新藥發(fā)現(xiàn)活動。在藥物發(fā)現(xiàn)過程的后期，MD模擬還可用于在進行更昂貴且更耗時的體外測試之前評估通過計算識別的小分子的質(zhì)量?？傊?#xff0c;MD 模擬在計算藥物設(shè)計中是一項寶貴工作。

在EFGR激酶上做分子動力學(xué)

表皮生長因子受體( EGFR ) 是與癌癥和炎癥相關(guān)的重要藥物靶點，制這種激酶是幾種已獲批準(zhǔn)的抗癌藥物的潛在作用機制，在本期中，我們使用這種激酶的 PDB 結(jié)構(gòu)3POZ（它與小分子抑制劑03P形成復(fù)合物）進行MD模擬。使用分子動力學(xué)引擎OpenMM進行分子動力學(xué)模擬。

首先準(zhǔn)備安裝包后導(dǎo)入依賴

conda?install?-q?-y?-c?conda-forge?mdtraj?openmm?openmmforcefields?openff-toolkit?pdbfixer?pypdb?rdkit
mamba?install?-q?-y?-c?conda-forge?openmmforcefieldsimport?copy
from?pathlib?import?Pathimport?requests
from?IPython.display?import?display
import?numpy?as?np
from?rdkit?import?Chem
from?rdkit.Chem?import?Draw
from?rdkit.Chem?import?AllChem
import?mdtraj?as?md
import?pdbfixer
import?openmm?as?mm
import?openmm.app?as?app
from?openmm?import?unit
from?openff.toolkit.topology?import?Molecule,?Topology
from?openmmforcefields.generators?import?GAFFTemplateGenerator

準(zhǔn)備文件地址，下載PDB文件

#?create?data?directory?if?not?exists
HERE?=?Path(_dh[-1])
DATA?=?HERE?/?"data"
DATA.mkdir(exist_ok=True)pdbid?=?"3POZ"
ligand_name?=?"03P"
pdb_path?=?DATA?/?f"{pdbid}.pdb"
pdb_url?=?f"https://files.rcsb.org/download/{pdbid}.pdb"r?=?requests.get(pdb_url)
r.raise_for_status()
with?open(pdb_path,?"wb")?as?f:f.write(r.content)

制備蛋白質(zhì)配體復(fù)合物

• 蛋白質(zhì)的制備：成功模擬的關(guān)鍵在于正確且完整的體系。從蛋白質(zhì)數(shù)據(jù)庫檢索到的晶體結(jié)構(gòu)通常會遺漏原子（主要是氫原子），并且可能包含非標(biāo)準(zhǔn)殘基。我們將使用 Python 中的PDBFixer來準(zhǔn)備蛋白質(zhì)結(jié)構(gòu)。然而， PDBFixer無法很好地處理共結(jié)晶配體，因此需要單獨準(zhǔn)備。

def?prepare_protein(pdb_file,?ignore_missing_residues=True,?ignore_terminal_missing_residues=True,?ph=7.0
):"""Use?pdbfixer?to?prepare?the?protein?from?a?PDB?file.?Hetero?atoms?such?as?ligands?areremoved?and?non-standard?residues?replaced.?Missing?atoms?to?existing?residues?are?added.Missing?residues?are?ignored?by?default,?but?can?be?included.Parameters----------pdb_file:?pathlib.Path?or?strPDB?file?containing?the?system?to?simulate.ignore_missing_residues:?bool,?optionalIf?missing?residues?should?be?ignored?or?built.ignore_terminal_missing_residues:?bool,?optionalIf?missing?residues?at?the?beginning?and?the?end?of?a?chain?should?be?ignored?or?built.ph:?float,?optionalpH?value?used?to?determine?protonation?state?of?residuesReturns-------fixer:?pdbfixer.pdbfixer.PDBFixerPrepared?protein?system."""fixer?=?pdbfixer.PDBFixer(str(pdb_file))fixer.removeHeterogens()??#?co-crystallized?ligands?are?unknown?to?PDBFixerfixer.findMissingResidues()??#?identify?missing?residues,?needed?for?identification?of?missing?atoms#?if?missing?terminal?residues?shall?be?ignored,?remove?them?from?the?dictionaryif?ignore_terminal_missing_residues:chains?=?list(fixer.topology.chains())keys?=?fixer.missingResidues.keys()for?key?in?list(keys):chain?=?chains[key[0]]if?key[1]?==?0?or?key[1]?==?len(list(chain.residues())):del?fixer.missingResidues[key]#?if?all?missing?residues?shall?be?ignored?ignored,?clear?the?dictionaryif?ignore_missing_residues:fixer.missingResidues?=?{}fixer.findNonstandardResidues()??#?find?non-standard?residuefixer.replaceNonstandardResidues()??#?replace?non-standard?residues?with?standard?onefixer.findMissingAtoms()??#?find?missing?heavy?atomsfixer.addMissingAtoms()??#?add?missing?atoms?and?residuesfixer.addMissingHydrogens(ph)??#?add?missing?hydrogensreturn?fixer#?prepare?protein?and?build?only?missing?non-terminal?residues
prepared_protein?=?prepare_protein(pdb_path,?ignore_missing_residues=False)

在這一步中，我們修復(fù)蛋白質(zhì)結(jié)構(gòu)中存在的問題（如缺失殘基、原子、非標(biāo)準(zhǔn)氨基酸、共結(jié)晶配體等），生成一個完整、標(biāo)準(zhǔn)化的拓?fù)淠Ｐ?#xff0c;用于分子動力學(xué)模擬。

• 準(zhǔn)備配體：蛋白質(zhì)制備完成后，我們將注意力轉(zhuǎn)向配體。同樣，我們需要添加氫原子，但也需要確保鍵序正確分配，因為某些 PDB 條目可能包含錯誤。我們使用 Python 包RDKit，這是一個開源化學(xué)信息學(xué)庫。我們將通過 SMILES 字符串向RDKit提供正確的質(zhì)子化狀態(tài)和鍵序。每個共結(jié)晶配體的不帶電異構(gòu)體 SMILES 字符串可以在其各自的PDB條目中找到。PDB 條目3POZ的配體名稱為03P。如果配體可能以其帶電形式或特定的互變異構(gòu)體形式結(jié)合，則需要將這些特性納入 SMILES 字符串中。

def?prepare_ligand(pdb_file,?resname,?smiles,?depict=True):"""Prepare?a?ligand?from?a?PDB?file?via?adding?hydrogens?and?assigning?bond?orders.?A?depictionof?the?ligand?before?and?after?preparation?is?rendered?in?2D?to?allow?an?inspection?of?theresults.?Huge?thanks?to?@j-wags?for?the?suggestion.Parameters----------pdb_file:?pathlib.PosixPathPDB?file?containing?the?ligand?of?interest.resname:?strThree?character?residue?name?of?the?ligand.smiles?:?strSMILES?string?of?the?ligand?informing?about?correct?protonation?and?bond?orders.depict:?bool,?optionalshow?a?2D?representation?of?the?ligandReturns-------prepared_ligand:?rdkit.Chem.rdchem.MolPrepared?ligand."""#?split?moleculerdkit_mol?=?Chem.MolFromPDBFile(str(pdb_file))rdkit_mol_split?=?Chem.rdmolops.SplitMolByPDBResidues(rdkit_mol)#?extract?the?ligand?and?remove?any?already?present?hydrogensligand?=?rdkit_mol_split[resname]ligand?=?Chem.RemoveHs(ligand)#?assign?bond?orders?from?templatereference_mol?=?Chem.MolFromSmiles(smiles)prepared_ligand?=?AllChem.AssignBondOrdersFromTemplate(reference_mol,?ligand)prepared_ligand.AddConformer(ligand.GetConformer(0))#?protonate?ligandprepared_ligand?=?Chem.rdmolops.AddHs(prepared_ligand,?addCoords=True)prepared_ligand?=?Chem.MolFromMolBlock(Chem.MolToMolBlock(prepared_ligand))#?2D?depictionif?depict:ligand_2d?=?copy.deepcopy(ligand)prepared_ligand_2d?=?copy.deepcopy(prepared_ligand)AllChem.Compute2DCoords(ligand_2d)AllChem.Compute2DCoords(prepared_ligand_2d)display(Draw.MolsToGridImage([ligand_2d,?prepared_ligand_2d],?molsPerRow=2,?legends=["original",?"prepared"]))#?return?ligandreturn?prepared_ligand

smiles?=?"CC(C)(O)CC(=O)NCCn1ccc2ncnc(Nc3ccc(Oc4cccc(c4)C(F)(F)F)c(Cl)c3)c12"
rdkit_ligand?=?prepare_ligand(pdb_path,?ligand_name,?smiles)

這一步是對我們即將加入其中的配體小分子，對照著原本蛋白中的配體的狀態(tài)進行清洗、結(jié)構(gòu)標(biāo)準(zhǔn)化和添加氫原子等處理，最終輸出一個可用于模擬的準(zhǔn)備好的 RDKit 分子對象。

• 合并蛋白質(zhì)和配體：下一步，我們要使用 Python 包MDTraj合并準(zhǔn)備好的蛋白質(zhì)和配體結(jié)構(gòu)。MDTraj可以處理準(zhǔn)備好的蛋白質(zhì)，它目前是一個PDBFixer分子，這種格式的拓?fù)浣Y(jié)構(gòu)和原子位置與OpenMM Modeller 的拓?fù)浣Y(jié)構(gòu)和位置相似，通常可以互換。然而，對于配體，我們需要進行一些轉(zhuǎn)換，因為它目前是一個RDKit分子。

def?rdkit_to_openmm(rdkit_mol,?name="LIG"):"""Convert?an?RDKit?molecule?to?an?OpenMM?molecule.Inspired?by?@hannahbrucemcdonald?and?@glass-w.Parameters----------rdkit_mol:?rdkit.Chem.rdchem.MolRDKit?molecule?to?convert.name:?strMolecule?name.Returns-------omm_molecule:?openmm.app.ModellerOpenMM?modeller?object?holding?the?molecule?of?interest."""#?convert?RDKit?to?OpenFFoff_mol?=?Molecule.from_rdkit(rdkit_mol)#?add?name?for?moleculeoff_mol.name?=?name#?add?names?for?atomselement_counter_dict?=?{}for?off_atom,?rdkit_atom?in?zip(off_mol.atoms,?rdkit_mol.GetAtoms()):element?=?rdkit_atom.GetSymbol()if?element?in?element_counter_dict.keys():element_counter_dict[element]?+=?1else:element_counter_dict[element]?=?1off_atom.name?=?element?+?str(element_counter_dict[element])#?convert?from?OpenFF?to?OpenMMoff_mol_topology?=?off_mol.to_topology()mol_topology?=?off_mol_topology.to_openmm()mol_positions?=?off_mol.conformers[0]#?convert?units?from??ngstr?m?to?nanometers#?since?OpenMM?works?in?nmmol_positions?=?mol_positions.to("nanometers")#?combine?topology?and?positions?in?modeller?objectomm_mol?=?app.Modeller(mol_topology,?mol_positions)return?omm_molomm_ligand?=?rdkit_to_openmm(rdkit_ligand,?ligand_name)

現(xiàn)在蛋白質(zhì)和配體都采用類似OpenMM 的格式，并且可以與MDTraj合并

def?merge_protein_and_ligand(protein,?ligand):"""Merge?two?OpenMM?objects.Parameters----------protein:?pdbfixer.pdbfixer.PDBFixerProtein?to?merge.ligand:?openmm.app.ModellerLigand?to?merge.Returns-------complex_topology:?openmm.app.topology.TopologyThe?merged?topology.complex_positions:?openmm.unit.quantity.QuantityThe?merged?positions."""#?combine?topologiesmd_protein_topology?=?md.Topology.from_openmm(protein.topology)??#?using?mdtraj?for?protein?topmd_ligand_topology?=?md.Topology.from_openmm(ligand.topology)??#?using?mdtraj?for?ligand?topmd_complex_topology?=?md_protein_topology.join(md_ligand_topology)??#?add?them?togethercomplex_topology?=?md_complex_topology.to_openmm()#?combine?positionstotal_atoms?=?len(protein.positions)?+?len(ligand.positions)#?create?an?array?for?storing?all?atom?positions?as?tupels?containing?a?value?and?a?unit#?called?OpenMM?Quantitiescomplex_positions?=?unit.Quantity(np.zeros([total_atoms,?3]),?unit=unit.nanometers)complex_positions[:?len(protein.positions)]?=?protein.positions??#?add?protein?positionscomplex_positions[len(protein.positions)?:]?=?ligand.positions??#?add?ligand?positionsreturn?complex_topology,?complex_positionscomplex_topology,?complex_positions?=?merge_protein_and_ligand(prepared_protein,?omm_ligand)

print("Complex?topology?has",?complex_topology.getNumAtoms(),?"atoms.")
#?NBVAL_CHECK_OUTPUT
Complex?topology?has?5561?atoms.

• MD模擬設(shè)置：

• 力場：常見的力場（例如 AMBER）包含氨基酸、核酸、水和離子的參數(shù)，通常提供多種選項，以滿足您的目標(biāo)。我們使用amber14-all.xmlOpenMM 自帶的力場文件，其中包含蛋白質(zhì)、DNA、RNA 和脂質(zhì)的參數(shù)。對于溶劑化，我們使用標(biāo)準(zhǔn)的三位點水模型 TIP3P。然而，配體的參數(shù)并未包含在內(nèi)。為了生成這些參數(shù)，我們可以使用通用AMBER力場( GAFF )，它由 Python 包OpenMM Forcefields實現(xiàn)。以下函數(shù)生成一個包含標(biāo)準(zhǔn)AMBER參數(shù)的力場對象，并根據(jù)需要額外包含小分子的參數(shù)。

def?generate_forcefield(rdkit_mol=None,?protein_ff="amber14-all.xml",?solvent_ff="amber14/tip3pfb.xml"
):"""Generate?an?OpenMM?Forcefield?object?and?register?a?small?molecule.Parameters----------rdkit_mol:?rdkit.Chem.rdchem.MolSmall?molecule?to?register?in?the?force?field.protein_ff:?stringName?of?the?force?field.solvent_ff:?stringName?of?the?solvent?force?field.Returns-------forcefield:?openmm.app.ForcefieldForcefield?with?registered?small?molecule."""forcefield?=?app.ForceField(protein_ff,?solvent_ff)if?rdkit_mol?is?not?None:gaff?=?GAFFTemplateGenerator(molecules=Molecule.from_rdkit(rdkit_mol,?allow_undefined_stereo=True))forcefield.registerTemplateGenerator(gaff.generator)return?forcefieldforcefield?=?generate_forcefield(rdkit_ligand)

• 系統(tǒng)：有了配置好的力場，我們現(xiàn)在可以使用OpenMM Modeller類來創(chuàng)建 MD 環(huán)境，這是一個包含復(fù)合物并充滿溶劑的模擬盒。標(biāo)準(zhǔn)溶劑是含有指定量離子的水。盒子的大小可以通過多種方式確定。我們用填充來定義它，這樣就會得到一個立方體盒子，其尺寸取決于復(fù)合物的最大尺寸。此步驟可能需要很長時間，大約幾分鐘，具體取決于您的硬件。

modeller?=?app.Modeller(complex_topology,?complex_positions)
modeller.addSolvent(forcefield,?padding=1.0?*?unit.nanometers,?ionicStrength=0.15?*?unit.molar)

有了溶劑化系統(tǒng)和力場，我們最終可以創(chuàng)建一個OpenMM 系統(tǒng)并設(shè)置模擬。除了系統(tǒng)之外，模擬還需要一個積分器。OpenMM積分器定義了一種通過積分運動方程來模擬系統(tǒng)的方法。所選的朗之萬積分器使用朗之萬方程。您可以在OpenMM 文檔中找到各種積分器的列表。

system?=?forcefield.createSystem(modeller.topology,?nonbondedMethod=app.PME)
integrator?=?mm.LangevinIntegrator(300?*?unit.kelvin,?1.0?/?unit.picoseconds,?2.0?*?unit.femtoseconds
)
simulation?=?app.Simulation(modeller.topology,?system,?integrator)
simulation.context.setPositions(modeller.positions)

• 執(zhí)行 MD 模擬：現(xiàn)在一切設(shè)置完畢，我們可以進行模擬了。我們需要設(shè)置起始位置并最小化系統(tǒng)能量，以獲得低能量的起始配置，這對于降低由于嚴(yán)重的原子碰撞導(dǎo)致模擬失敗的可能性至關(guān)重要。能量最小化的系統(tǒng)已保存。

simulation.minimizeEnergy()
with?open(DATA?/?"topology.pdb",?"w")?as?pdb_file:app.PDBFile.writeFile(simulation.topology,simulation.context.getState(getPositions=True,?enforcePeriodicBox=True).getPositions(),file=pdb_file,keepIds=True,)

最小化完成后，我們就可以進行分子動力學(xué)模擬了。在本講座中，我們將進行一個簡短的模擬來進行說明。用于研究的模擬通?？缭綆准{秒，甚至幾微秒。我們將僅模擬100皮秒的分子動力學(xué)，對應(yīng)5萬步，每步2飛秒。我們每10皮秒（5000步）保存一次分子“快照”，總共10幀。結(jié)果保存在一個.xtc文件中，該文件包含給定時間點所有原子的坐標(biāo)。結(jié)合之前編寫的能量最小化體系的PDB文件，它為我們提供了后續(xù)分析所需的所有信息。如果您有性能良好的 GPU，也可以增加模擬時間。

#?output?settings
if?on_colab:steps?=?50000??#?corresponds?to?100?pswrite_interval?=?5000??#?write?every?10?pslog_interval?=?2500??#?log?progress?to?stdout?every?5?ps
else:steps?=?10??#?corresponds?to?20?fswrite_interval?=?1??#?write?every?2?fslog_interval?=?1??#?log?progress?to?stdout?every?2?fs
simulation.reporters.append(md.reporters.XTCReporter(file=str(DATA?/?"trajectory.xtc"),?reportInterval=write_interval)
)
simulation.reporters.append(app.StateDataReporter(sys.stdout,log_interval,step=True,potentialEnergy=True,temperature=True,progress=True,remainingTime=True,speed=True,totalSteps=steps,separator="\t",)
)

系統(tǒng)中所有粒子的速度都是根據(jù)給定溫度下的分布隨機選擇的。我們選擇了 300 開爾文，比室溫高出幾度。隨機種子會生成，但也可以明確指定以重現(xiàn)結(jié)果。

然后按照之前定義的步驟進行模擬。

simulation.context.setVelocitiesToTemperature(300?*?unit.kelvin)
simulation.step(steps)??#?perform?the?simulation
#"Progress?(%)"?"Step"??"Potential?Energy?(kJ/mole)"????"Temperature?(K)"???????"Speed?(ns/day)"????????"Time?Remaining"
10.0%???1???????-992629.2532507611??????280.59188983053207??????0???????--
20.0%???2???????-978514.3140385933??????250.1140053691449???????1.88????0:00
30.0%???3???????-961427.9816461714??????214.09328454071326??????1.92????0:00
40.0%???4???????-944609.5657245225??????180.29949178380076??????1.9?????0:00
50.0%???5???????-938426.3768274506??????165.17778679274846??????1.89????0:00
60.0%???6???????-930096.1319021038??????149.4126629808844???????1.89????0:00
70.0%???7???????-926924.8985479049??????143.336941762251????????1.9?????0:00
80.0%???8???????-930284.4123182578??????149.77757492192208??????1.91????0:00
90.0%???9???????-932318.9156875424??????155.78884925604478??????1.9?????0:00
100.0%??10??????-937770.7580560565??????167.02287919129623??????1.9?????0:00

讓我們來看看可視化結(jié)果 第一步，我們加載拓?fù)?( .pdb) 和軌跡 ( .xtc) 文件，以創(chuàng)建表示模擬系統(tǒng)的 MDAnalysis universe。NGL Viewer提供預(yù)構(gòu)建功能來可視化加載在常用庫（如 MD Analysis）中的模擬。

from?pathlib?import?Pathimport?pandas?as?pd
import?numpy?as?np
import?matplotlib.pyplot?as?plt
import?nglview?as?nv
import?MDAnalysis?as?mda
from?MDAnalysis.analysis?import?rms,?diffusionmap,?align
from?MDAnalysis.analysis.distances?import?dist
from?MDAnalysis.analysis.hydrogenbonds.hbond_analysis?import?HydrogenBondAnalysis?as?HBAmd_universe?=?mda.Universe(str(DATA?/?"topology.pdb"),?str(DATA?/?"trajectory.xtc"))
view?=?nv.show_mdanalysis(md_universe)
view
view.render_image(trim=True,?factor=2,?transparent=True);
view._display_image()

這一步可以看見動態(tài)的結(jié)果哦

• 隨時間變化的 RMSD 蛋白質(zhì)和配體的 RMSD? 動畫中蛋白質(zhì)和配體清晰可辨。然而，很難評估配體或蛋白質(zhì)是否真的顯著改變了它們的位置/構(gòu)象。因此，我們將使用更完善的方法——計算均方根差 (RMSD)——來評估相對運動。因此，我們將重點關(guān)注系統(tǒng)的不同部分（骨架、蛋白質(zhì)和配體）。配體最好通過其獨特的殘基名稱來選擇。在本例中，我們知道它是抑制劑03P，但我們也可以通過檢查拓?fù)湮募驅(qū)⒐鈽?biāo)懸停在上面動畫中的配體原子上來找到它的名稱。為了繪制 RMSD，我們編寫了一個函數(shù)，對從宇宙中選定的原子組進行 RMSD 計算并返回一個pandas Dataframe。

ligand_name?=?"03P"def?rmsd_for_atomgroups(universe,?selection1,?selection2=None):"""Calulate?the?RMSD?for?selected?atom?groups.Parameters----------universe:?MDAnalysis.core.universe.UniverseMDAnalysis?universe.selection1:?strSelection?string?for?main?atom?group,?also?used?during?alignment.selection2:?list?of?str,?optionalSelection?strings?for?additional?atom?groups.Returns-------rmsd_df:?pandas.core.frame.DataFrameDataFrame?containing?RMSD?of?the?selected?atom?groups?over?time."""universe.trajectory[0]ref?=?universermsd_analysis?=?rms.RMSD(universe,?ref,?select=selection1,?groupselections=selection2)rmsd_analysis.run()columns?=?[selection1,?*selection2]?if?selection2?else?[selection1]rmsd_df?=?pd.DataFrame(np.round(rmsd_analysis.results.rmsd[:,?2:],?2),?columns=columns)rmsd_df.index.name?=?"frame"return?rmsd_df

運行該函數(shù)來計算模擬過程中選定分子的 RMSD：全蛋白質(zhì)、蛋白質(zhì)骨架原子、配體：

在計算過程中，分子將沿著蛋白質(zhì)的主鏈疊加到第一幀

rmsd?=?rmsd_for_atomgroups(md_universe,?"backbone",?["protein",?f"resname?{ligand_name}"])
rmsd.head()

然后，我們繪制結(jié)果。

rmsd.plot(title="RMSD?of?protein?and?ligand")
plt.ylabel("RMSD?(?)");

從中我們能得到什么結(jié)果：

1. 在模擬過程中，主鏈僅發(fā)生了輕微的構(gòu)象變化，這并不奇怪，因為蛋白質(zhì)是相對于主鏈疊加的。

2. 與僅計算蛋白質(zhì)主鏈相比，將蛋白質(zhì)側(cè)鏈納入 RMSD 計算中會導(dǎo)致更多的運動。

3. 在我們的模擬中，配體的運動范圍與蛋白質(zhì)側(cè)鏈相同。然而，情況并非總是如此，因為配體可以更自由地移動，而且除了構(gòu)象變化之外，還可能涉及定向運動。如果配體的 RMSD 值非常高，則很可能出現(xiàn)這種情況。

計算所有幀（成對）之間的差異，并在熱圖中將其可視化。

def?RMSD_dist_frames(universe,?selection):"""Calculate?the?RMSD?between?all?frames?in?a?matrix.Parameters----------universe:?MDAnalysis.core.universe.UniverseMDAnalysis?universe.selection:?strSelection?string?for?the?atomgroup?to?be?investigated,?also?used?during?alignment.Returns-------array:?np.ndarrayNumpy?array?of?RMSD?values."""pairwise_rmsd?=?diffusionmap.DistanceMatrix(universe,?select=selection)pairwise_rmsd.run()return?pairwise_rmsd.results.dist_matrixdist_matrix_protein?=?RMSD_dist_frames(md_universe,?"protein")
dist_matrix_ligand?=?RMSD_dist_frames(md_universe,?f"resname?{ligand_name}")#?get?maximal?distance?between?all?frames
max_dist?=?max(np.amax(dist_matrix_ligand),?np.amax(dist_matrix_protein))fig,?ax?=?plt.subplots(1,?2)
fig.suptitle("RMSD?between?the?frames")#?protein?image
img1?=?ax[0].imshow(dist_matrix_protein,?cmap="viridis",?vmin=0,?vmax=max_dist)
ax[0].title.set_text("protein")
ax[0].set_xlabel("frames")
ax[0].set_ylabel("frames")#?ligand?image
img2?=?ax[1].imshow(dist_matrix_ligand,?cmap="viridis",?vmin=0,?vmax=max_dist)
ax[1].title.set_text("Ligand")
ax[1].set_xlabel("frames")fig.colorbar(img1,?ax=ax,?orientation="horizontal",?fraction=0.1,?label="RMSD?(?)");

從熱圖可以看出，蛋白質(zhì)在不同狀態(tài)之間持續(xù)均勻地移動。配體在不同狀態(tài)之間跳躍。有時它可以反復(fù)達(dá)到相似的構(gòu)象。