{ "cells": [ { "cell_type": "markdown", "id": "c4185e43-9c43-400c-9688-7b3831dacecb", "metadata": {}, "source": [ "## Sphere Model: Brownian Dynamics " ] }, { "cell_type": "markdown", "id": "bc7ad61e-e7f4-4b91-8df2-21995cd485ed", "metadata": {}, "source": [ "### Import the necessary packages" ] }, { "cell_type": "code", "execution_count": 1, "id": "d7dfe668-ce05-48b3-98ba-a5b73e0df96d", "metadata": {}, "outputs": [], "source": [ "import parmed as pmd\n", "import json\n", "import sys\n", "from sys import platform\n", "import mdtraj as md\n", "import simtk.openmm.app as mmapp\n", "import mdtraj.reporters\n", "import simtk.unit as u\n", "import random\n", "from openNucleome import OpenNucleome\n", "from openNucleome.utils import coor_transformation, final_frame\n", "import warnings\n", "import numpy as np" ] }, { "cell_type": "markdown", "id": "d8c3b70a-21f3-4228-bb57-945965ebee06", "metadata": {}, "source": [ "### Important parameters\n", "We defined the important parameters in the next block. First, we set the transition probability between dP particles and P particles as 0.2, and the transition frequency as 4000. When creating the system, we set type 6, 7 as the dP and P particles, respectively, so we kept using 6, 7 here. In this example, we ran a simulation with total length of 100,000 steps, and output one configuration and the energy every 2000 steps." ] }, { "cell_type": "code", "execution_count": 2, "id": "5a6e5ed3-5b95-4568-af71-336b5564cafa", "metadata": {}, "outputs": [], "source": [ "prob_P_dP = 0.2 # Transition probability from P to dP\n", "prob_dP_P = 0.2 # Transition probability from dP to P\n", "transition_freq = 4000\n", "sampling_freq = 2000\n", "dP_type = 6\n", "P_type = 7\n", "total_steps = 100000" ] }, { "cell_type": "markdown", "id": "95f76ca2-43ed-4768-9fce-257348ea81c0", "metadata": {}, "source": [ "### Initialize the system (A larger damping coefficient for the Brownian Dynamics)\n", "We first set up an example \"model\" of class \"OpenNucleome\" with the conserved temperature, damping coefficient, timestep, and the mass scale. In this folder, we also included the initial configuration \"human.pdb\" used for the simulation and created a system according to the initial configuration.\n", "\n", "In this example, we freezed all the lamina beads, and would not consider the dynamics of the membrane, and that is why we set \"False\" (off) for membrane dynamics. Consequently, there was also no need to set the bond between lamina beads, so we set \"None\" for the variable \"lam_bond\"." ] }, { "cell_type": "code", "execution_count": 3, "id": "b0b6fb41-d073-4be8-8832-24896e0f068f", "metadata": {}, "outputs": [], "source": [ "model = OpenNucleome(1.0, 100, 0.005, 1.0) # 1.0: temperature (LJ reduced unit); \n", " # 0.1: damping coefficient (LJ reduced unit);\n", " # 0.005: timestep (LJ reduced unit);\n", " # 1.0: mass_scale\n", " \n", "PDB_file = \"initial_configs/human_1.pdb\" #The initial configuration\n", "\n", "# Generate new elements and construct topology as well\n", "# flag_membrane: True for including lamina dynamics, False for excluding lamina dynamics;\n", "# lam_bond: A file contains the lamina bond when membrane_dynamics is on.\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "model.create_system(PDB_file, flag_membrane = False, lam_bond = None) " ] }, { "cell_type": "markdown", "id": "3eee92a7-e149-45ef-b201-7b02fa5948c6", "metadata": {}, "source": [ "### Add the force field\n", "In this example, we loaded default settings for the interactions between chromosomes and chromosomes, and chromosomes and nuclear landmarks. All the types of potential can be found in \"Section: Energy function of the whole nucleus model\" in Supporting Information. According to the order of added potential, the index of speckle-speckle potential is 6, and this index is needed when transiting the speckle particle between type 6 and type 7.\n", "\n", "Because during the simulations, we would transit the dP speckle and P speckle particles, here, we logged the start index and end index of speckle particle, and computed the number of speckle particles." ] }, { "cell_type": "code", "execution_count": 4, "id": "9c8cf2a5-8071-4e29-82d5-c087905883bf", "metadata": {}, "outputs": [], "source": [ "# Add the default force field\n", "model.load_default_settings()\n", "\n", "index_spec_spec_potential = 6 #The 6-th potential is speckle-speckle interaction in the default settings\n", "start_spec_index = model.N_chr_nuc+1\n", "end_spec_index = model.N_chr_nuc_spec+1\n", "N_spec = end_spec_index-start_spec_index #Currently, the number of speckle is 1600." ] }, { "cell_type": "markdown", "id": "d5d920a1-35a3-4258-a774-2d351b3efc46", "metadata": {}, "source": [ "### Perform the simulation\n", "We first created the simulation with a specific Platform, and here, we used \"CUDA\" but users can also use \"CPU\", \"Reference\", and \"OpenCL\" according to their hardware. Before the simulation, we minimized the energy to make the system much more reasonable and stable. After randomly setting velocity, we started our simulation with a total length of 100,000 steps and output the configuration and energy every 2000 steps, and change the speckle types every 4000 steps as we mentioned in the previous blocks." ] }, { "cell_type": "code", "execution_count": 5, "id": "f64c55d0-1db4-4897-a4be-1a774f26281e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "#\"Progress (%)\"\t\"Step\"\t\"Potential Energy (kJ/mole)\"\t\"Kinetic Energy (kJ/mole)\"\t\"Temperature (K)\"\t\"Time Remaining\"\n", "2.0%\t2000\t99669.37249189641\t93484.31840286203\t119.85097090889614\t--\n", "4.0%\t4000\t102879.7194162421\t93650.98329636063\t120.06464256681159\t12:11\n", "6.0%\t6000\t104829.51531281356\t93560.6536394227\t119.94883600940348\t11:58\n", "8.0%\t8000\t105119.10886859841\t93785.34285541113\t120.23689737786663\t11:42\n", "10.0%\t10000\t105999.76281040246\t93984.31355069151\t120.49198648178462\t11:27\n", "12.0%\t12000\t105999.28209041196\t93586.10689548533\t119.98146819311104\t11:11\n", "14.0%\t14000\t107031.57428997077\t93881.41950974891\t120.36007183643666\t10:56\n", "16.0%\t16000\t106615.33266635888\t93369.33133546957\t119.7035524765387\t10:41\n", "18.0%\t18000\t106557.47937161615\t94058.32191220552\t120.58686842706025\t10:26\n", "20.0%\t20000\t107007.09729298094\t93408.50286975966\t119.75377209087473\t10:10\n", "22.0%\t22000\t107926.75086174044\t94138.0784156327\t120.6891197408138\t9:55\n", "24.0%\t24000\t107391.15834572143\t93151.26510476637\t119.42398206388837\t9:40\n", "26.0%\t26000\t107501.08014869876\t94274.11024334499\t120.8635184731827\t9:24\n", "28.0%\t28000\t107569.04196787113\t94384.1702547747\t121.00462020503947\t9:09\n", "30.0%\t30000\t108045.97777092032\t94601.02906711765\t121.282642654724\t8:54\n", "32.0%\t32000\t107432.52383930977\t93589.96460811926\t119.9864139488545\t8:39\n", "34.0%\t34000\t108329.11488470988\t93621.14913681203\t120.02639387388383\t8:23\n", "36.0%\t36000\t107958.7920328365\t93807.04607312058\t120.26472184895187\t8:08\n", "38.0%\t38000\t108615.67340018699\t93614.86099439148\t120.01833220121839\t7:53\n", "40.0%\t40000\t108189.84506806819\t93740.44910207292\t120.17934162929558\t7:37\n", "42.0%\t42000\t107939.65143493259\t93834.05573605302\t120.29934941410589\t7:22\n", "44.0%\t44000\t107817.74457492097\t94328.42213398592\t120.93314868428571\t7:07\n", "46.0%\t46000\t107645.46570848562\t93837.58325264072\t120.30387184412577\t6:52\n", "48.0%\t48000\t108009.67809897116\t93728.81099086833\t120.16442106345902\t6:36\n", "50.0%\t50000\t108136.9096828193\t93568.82887522175\t119.9593170180332\t6:21\n", "52.0%\t52000\t107895.9760849019\t93795.82888007448\t120.25034091854111\t6:06\n", "54.0%\t54000\t107886.52785448956\t94409.24909202677\t121.03677236751204\t5:51\n", "56.0%\t56000\t108117.69439003256\t93699.94336879319\t120.12741151370089\t5:35\n", "58.0%\t58000\t108090.49330945831\t93610.20515342659\t120.0123632101665\t5:20\n", "60.0%\t60000\t107685.93326105399\t94339.93374000705\t120.94790706496308\t5:05\n", "62.0%\t62000\t108413.78749640704\t93574.87401785952\t119.96706715438546\t4:50\n", "64.0%\t64000\t108240.0185759508\t94063.73526929892\t120.59380858680689\t4:34\n", "66.0%\t66000\t108650.74974339592\t93739.3187409187\t120.17789245702663\t4:19\n", "68.0%\t68000\t108106.430386827\t93157.33329748194\t119.43176175144505\t4:04\n", "70.0%\t70000\t108253.03044933111\t94041.39083856613\t120.56516205267707\t3:48\n", "72.0%\t72000\t107594.28937622189\t93723.37105782727\t120.15744683218914\t3:33\n", "74.0%\t74000\t108512.84246920235\t93776.25733124021\t120.22524933986725\t3:18\n", "76.0%\t76000\t107968.59366912623\t93890.56791089648\t120.37180048546071\t3:03\n", "78.0%\t78000\t108024.27839576988\t93862.2149239762\t120.33545072040258\t2:47\n", "80.0%\t80000\t107697.96213950164\t93773.7042630089\t120.22197619511371\t2:32\n", "82.0%\t82000\t108108.68537563692\t94185.73338566914\t120.75021548954557\t2:17\n", "84.0%\t84000\t108250.88620647619\t93754.43863044186\t120.19727681442845\t2:02\n", "86.0%\t86000\t108259.21231578765\t93841.62919263162\t120.30905891554517\t1:46\n", "88.0%\t88000\t108737.65904237027\t93828.09609682369\t120.29170889685012\t1:31\n", "90.0%\t90000\t108590.35111284332\t94000.81867313881\t120.51314676817357\t1:16\n", "92.0%\t92000\t108451.29548567365\t93910.97405537797\t120.39796205213196\t1:01\n", "94.0%\t94000\t108650.86069591212\t94028.24587136641\t120.54830963177517\t0:45\n", "96.0%\t96000\t108038.8927258845\t94171.9106591098\t120.73249415160302\t0:30\n", "98.0%\t98000\t108029.91821528843\t93447.38024299666\t119.80361457791906\t0:15\n", "100.0%\t100000\t108042.33377743122\t94088.84286342005\t120.62599761683903\t0:00\n" ] } ], "source": [ "#model.save_system(\"model_before_simulation_0.xml\")\n", "\n", "simulation = model.create_simulation(platform_type = \"CUDA\") # Users can also use CPU, Reference, OpenCL.\n", "simulation.context.setPositions(model.chr_positions)\n", "\n", "simulation.minimizeEnergy()\n", "\n", "simulation.reporters.append(mdtraj.reporters.DCDReporter('results/step_100000.dcd', sampling_freq))\n", "\n", "def setVelocity(context):\n", " sigma = u.sqrt(1.0*u.kilojoule_per_mole / model.chr_system.getParticleMass(1))\n", " velocs = u.Quantity(1.0 * np.random.normal(size=(model.chr_system.getNumParticles(), 3)),\n", " u.meter) * (sigma / u.meter)\n", " context.setVelocities(velocs)\n", "setVelocity(simulation.context)\n", "\n", "simulation.reporters.append(mmapp.statedatareporter.StateDataReporter(sys.stdout, sampling_freq, step=True,\n", " potentialEnergy=True, kineticEnergy=True, temperature=True, progress=True,\n", " remainingTime=True, separator='\\t', totalSteps = total_steps))\n", "\n", "for i in range(total_steps//transition_freq):\n", " simulation.step(transition_freq)\n", " # Change the type of speckles every 4000 steps, non-equilibrium scheme.\n", "\n", " # Do the chemical modification, and change the spec-spec potential on the fly.\n", " for j in np.random.randint(start_spec_index-1, end_spec_index-1, N_spec): \n", "\n", " if model.compart_type[j] == dP_type-1:\n", " model.compart_type[j] = P_type-1 if random.random() < prob_dP_P else dP_type-1\n", " else:\n", " model.compart_type[j] = dP_type-1 if random.random() < prob_P_dP else P_type-1\n", "\n", " # Update the context after changing the type of speckles.\n", " for m in range(model.chr_system.getNumParticles()):\n", " model.chr_system.getForce(index_spec_spec_potential).setParticleParameters(m, [model.compart_type[m]])\n", " model.chr_system.getForce(index_spec_spec_potential).updateParametersInContext(simulation.context)\n", "\n", "# Keep the final result of bead types in case constructing the configuration for the continuous simulation.\n", "np.savetxt('results/compt_final_frame.txt', (np.array(model.compart_type)+1).reshape((-1,1)), fmt='%d')" ] }, { "cell_type": "markdown", "id": "ef20c542-a6d7-4e35-b398-2cc17e3ee149", "metadata": {}, "source": [ "### Post-process the output files\n", "\n", "We converted the trajectory with the default OpenMM unit to the nex trajectory with the reduced unit, and saved the final configuration according to the last frame." ] }, { "cell_type": "code", "execution_count": 6, "id": "be56b164-2e75-4c03-9636-22e023298fb5", "metadata": {}, "outputs": [], "source": [ "coor_transformation('results/step_100000.dcd', 'results/reduced_traj.dcd')\n", "final_frame('results/step_100000.dcd', 'results/compt_final_frame.txt', 'results/final_config.pdb')" ] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:anaconda3-openmm]", "language": "python", "name": "conda-env-anaconda3-openmm-py" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.15" } }, "nbformat": 4, "nbformat_minor": 5 }