{ "cells": [ { "cell_type": "markdown", "id": "c4185e43-9c43-400c-9688-7b3831dacecb", "metadata": {}, "source": [ "## Nuclear Deformation Simulation" ] }, { "cell_type": "markdown", "id": "41a7c311-1202-428e-b9ce-ed77aa648d0e", "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\n", "import pandas as pd" ] }, { "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\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 moved all the lamina beads, and would consider the dynamics of the membrane, and that is why we set \"True\" (on) for membrane dynamics. Consequently, we need to designate the bonds between specific lamina beads, so we included a text file, which logs all the pairs linked by the bonds, for the variable \"membrane_bond\"." ] }, { "cell_type": "code", "execution_count": 3, "id": "b0b6fb41-d073-4be8-8832-24896e0f068f", "metadata": {}, "outputs": [], "source": [ "model = OpenNucleome(1.0, 0.1, 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 = \"strength_one/human.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 = True, lam_bond = 'input_params/lamina_bond.txt') " ] }, { "cell_type": "markdown", "id": "3eee92a7-e149-45ef-b201-7b02fa5948c6", "metadata": {}, "source": [ "### Add the force field\n", "Different from the sphere nucleus situation, we used a csv file to log all the flags of interactions between chromosomes and chromosomes, or chromosomes and nuclear landmarks.\n", "\n", "The column from \"bond\" to \"lam_squeeze\" is boolen, which means the switches of the corresponding interactions. For example, if the spec-chrom is True, it means in the simulation, the interaction between speckle and chromosome will be present; if the inter is False, it means we exclude the interchromosomal interactions between chromosome. Remember some potentials require the corresponding parameter files, such as ideal requires ideal_param_file and spec-chrom requires the chr_spec_param, so if you turn those potentials on, you should also include the corresponding parameter files.\n", "\n", "In this example, we turn every interaction on, and set the strength of force squeezing the nucleus as 1.0 (k) when loading the customized force field. The users can set their specific strengths.\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": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
bondanglesoftcoreidealcomptinterspec-specspec-chromnuc-nucnuc-spec...lam-chromhard-walllam-lamlam_squeezeideal_param_filecompt_param_fileinterchr_param_filechr_spec_paramchr_nuc_paramchr_lam_param
chromosomeTrueTrueTrueTrueTrueTrueNaNNaNNaNNaN...NaNNaNNaNNaNideal_param_file.txtcompt_param_file.txtinterchr_param_file.txtNaNNaNNaN
speckleNaNNaNNaNNaNNaNNaNTrueTrueNaNNaN...NaNNaNNaNNaNNaNNaNNaNchr_spec_param.txtNaNNaN
nucleolusNaNNaNNaNNaNNaNNaNNaNNaNTrueTrue...NaNNaNNaNNaNNaNNaNNaNNaNchr_nuc_param.txtNaN
laminaNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN...TrueTrueTrueTrueNaNNaNNaNNaNNaNchr_lam_param.txt
\n", "

4 rows × 21 columns

\n", "
" ], "text/plain": [ " bond angle softcore ideal compt inter spec-spec spec-chrom \\\n", "chromosome True True True True True True NaN NaN \n", "speckle NaN NaN NaN NaN NaN NaN True True \n", "nucleolus NaN NaN NaN NaN NaN NaN NaN NaN \n", "lamina NaN NaN NaN NaN NaN NaN NaN NaN \n", "\n", " nuc-nuc nuc-spec ... lam-chrom hard-wall lam-lam lam_squeeze \\\n", "chromosome NaN NaN ... NaN NaN NaN NaN \n", "speckle NaN NaN ... NaN NaN NaN NaN \n", "nucleolus True True ... NaN NaN NaN NaN \n", "lamina NaN NaN ... True True True True \n", "\n", " ideal_param_file compt_param_file \\\n", "chromosome ideal_param_file.txt compt_param_file.txt \n", "speckle NaN NaN \n", "nucleolus NaN NaN \n", "lamina NaN NaN \n", "\n", " interchr_param_file chr_spec_param chr_nuc_param \\\n", "chromosome interchr_param_file.txt NaN NaN \n", "speckle NaN chr_spec_param.txt NaN \n", "nucleolus NaN NaN chr_nuc_param.txt \n", "lamina NaN NaN NaN \n", "\n", " chr_lam_param \n", "chromosome NaN \n", "speckle NaN \n", "nucleolus NaN \n", "lamina chr_lam_param.txt \n", "\n", "[4 rows x 21 columns]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Add the default force field\n", "# Add the customized force field\n", "# In this customized setting, I logged all the flags for different interactions\n", "# (True means on, False means off, and NaN means N/A)\n", "# and all the corresponding files required for interactions. \n", "force_field = pd.read_csv('input_params/input.csv', sep=' ', header=0, index_col=0)\n", "model.load_customized_settings(force_field, 'input_params', k = 1.0)\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\n", "\n", "force_field" ] }, { "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\t123136.44445905171\t100752.69234185581\t114.52055593216552\t--\n", "4.0%\t4000\t138145.93252829043\t99854.70038122084\t113.4998532971868\t9:06\n", "6.0%\t6000\t142575.6182367836\t104465.90880812544\t118.741183730089\t8:38\n", "8.0%\t8000\t141094.42452352808\t107327.50563145665\t121.99381799170929\t8:21\n", "10.0%\t10000\t139697.16552712512\t107342.47070395519\t122.01082804258003\t8:13\n", "12.0%\t12000\t137654.04464475974\t106994.70033531352\t121.61553482482218\t8:03\n", "14.0%\t14000\t138386.55393252496\t105782.62134612219\t120.2378251433289\t7:53\n", "16.0%\t16000\t137409.19082230702\t106199.98460544148\t120.71222112592687\t7:42\n", "18.0%\t18000\t138944.10917304258\t106093.39691125305\t120.59106821467195\t7:32\n", "20.0%\t20000\t137593.22835381792\t106262.0393691768\t120.78275567816587\t7:21\n", "22.0%\t22000\t138695.92867257126\t105989.21260070454\t120.47264710962654\t7:10\n", "24.0%\t24000\t138082.96422167047\t105636.71725733916\t120.07198324896527\t6:59\n", "26.0%\t26000\t138577.68014618923\t105942.89960543014\t120.42000543978695\t6:48\n", "28.0%\t28000\t138547.97569732778\t105614.2894846135\t120.04649071927538\t6:37\n", "30.0%\t30000\t139048.79707005448\t105953.13915726854\t120.43164422721559\t6:26\n", "32.0%\t32000\t139369.4478256774\t106096.93739319986\t120.5950925037643\t6:15\n", "34.0%\t34000\t139349.70734011347\t105551.45826302844\t119.9750736061584\t6:04\n", "36.0%\t36000\t138652.69678783504\t106307.39016938012\t120.83430366889978\t5:53\n", "38.0%\t38000\t138994.0196299976\t105924.91785091767\t120.39956647706012\t5:42\n", "40.0%\t40000\t138711.95389602106\t106045.77361075555\t120.53693718629414\t5:31\n", "42.0%\t42000\t139755.55680584465\t106157.32762501454\t120.66373506566056\t5:20\n", "44.0%\t44000\t139418.1188406263\t105861.9173488944\t120.32795695131931\t5:09\n", "46.0%\t46000\t139318.77307771594\t106084.84309645879\t120.58134552039144\t4:58\n", "48.0%\t48000\t139156.1582248191\t106000.44885852029\t120.48541880298527\t4:47\n", "50.0%\t50000\t139394.24871664512\t106154.93144347834\t120.66101144572286\t4:36\n", "52.0%\t52000\t138976.02488041663\t105917.33440921981\t120.39094675743841\t4:24\n", "54.0%\t54000\t139636.7751483859\t105824.53815483821\t120.2854698873599\t4:13\n", "56.0%\t56000\t139109.49119912702\t105962.1255783129\t120.44185864342542\t4:02\n", "58.0%\t58000\t139591.83060822415\t105548.44699685683\t119.9716508502199\t3:51\n", "60.0%\t60000\t139158.00653002597\t105773.10164969457\t120.22700457961139\t3:40\n", "62.0%\t62000\t139141.23690284294\t106340.19938493858\t120.87159626642863\t3:29\n", "64.0%\t64000\t138657.07116640714\t106162.69666414094\t120.66983778441957\t3:18\n", "66.0%\t66000\t138872.56066384647\t105709.55038510424\t120.154769029554\t3:07\n", "68.0%\t68000\t139552.46901895298\t106237.45014421445\t120.75480633361563\t2:56\n", "70.0%\t70000\t140314.6318439555\t105998.06802734527\t120.48271263103683\t2:45\n", "72.0%\t72000\t139878.4964976344\t105824.37540777247\t120.28528490089441\t2:34\n", "74.0%\t74000\t140049.54976571468\t105978.09999204957\t120.46001596209643\t2:23\n", "76.0%\t76000\t139040.28804717027\t105698.80343463195\t120.14255350742695\t2:12\n", "78.0%\t78000\t139441.61604379857\t105796.50827126215\t120.25360971791694\t2:01\n", "80.0%\t80000\t139103.81283076043\t105733.9975154245\t120.1825568622178\t1:50\n", "82.0%\t82000\t139735.49375639576\t105520.89870734174\t119.94033808470546\t1:39\n", "84.0%\t84000\t139625.44875567633\t105857.33796241746\t120.32275179130842\t1:28\n", "86.0%\t86000\t139302.36985220472\t106414.25964148348\t120.95577686304458\t1:17\n", "88.0%\t88000\t139197.9629622838\t105156.448308775\t119.52608550957876\t1:06\n", "90.0%\t90000\t138969.3590577207\t105676.13902677708\t120.1167920063982\t0:55\n", "92.0%\t92000\t138261.92996336846\t105622.26065619843\t120.05555115210468\t0:44\n", "94.0%\t94000\t138923.98136101887\t106060.5908358773\t120.55377918643204\t0:33\n", "96.0%\t96000\t138863.37451388055\t105655.18511185612\t120.0929747373017\t0:22\n", "98.0%\t98000\t138649.21396001056\t106102.68049923926\t120.60162040571038\t0:11\n", "100.0%\t100000\t138531.91066385998\t106315.62779116184\t120.84366696236692\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": "d1d36075-fdef-4ef0-95a9-d4d0bc3a6b82", "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": "6a73e27e-6326-4e45-9da3-299907a60fa7", "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 }