Course: BT2042: Fundamentals of Biophysical Chemistry
Authors: Ritvik Pandey (BS21B028), S Deepak Kumar (BS21B013), Rahul Srinivas (BS21B029)
This repository contains the simulation code and data analysis for the unfolding of a 2D 16-mer polymer. The simulations start from a fully folded state and explore the unfolding dynamics under various non-covalent interaction energy values (-0.25, -0.75, -1.25, and -2.5 kT).
The simulation utilizes a Monte Carlo approach to sample the conformational space of the polymer.
- Residue Mapping: The protein structure is tracked using a 2D NumPy array storing the coordinates of each of the 16 residues.
-
Move Sets: * End Moves: Available only for residues 1 and 16. Allows movement to 4 adjacent diagonal coordinates if no overlap occurs and bond length is conserved.
- Corner Moves: Allows intermediate residues to move diagonally if the previous and next residues are fixed, maintaining structural integrity.
-
Metropolis Criterion: Move acceptance is dictated by the energy difference between states (
$E_{i+1} - E_i$ ). Favorable moves are accepted; unfavorable moves are accepted with a probability of$e^{(-\Delta E/kT)}$ . - Simulation Length: 100,000 steps per run, averaged over 50 individual molecules for each energy state.
The script tracks multiple structural and thermodynamic parameters to monitor unfolding:
-
Radius of Gyration (
$R_g$ ): Measures the compactness of the polymer configuration around its center of mass. - End-to-End Distance: The spatial distance between the 1st and 16th residue.
-
Interaction Energy: Calculated as
$E = N imes E_S$ , where$N$ is the number of native non-covalent interactions. -
Microstates: Sampled and identified uniquely based on a combined score of Interaction Energy,
$R_g$ , and End-to-End distance. -
Free Energy (
$\Delta G$ ): Evaluated using$\Delta G = \Delta U - T\Delta S$ , utilizing the number of sampled microstates to approximate entropy.
- Energy Dependence: As interaction energies become more negative (stronger), the polymer requires significantly more turns to unfold. For example, at -0.25 kT the mean unfolding turns is ~279, whereas at -2.50 kT it jumps to ~995.
-
Free Energy Profiles: The reaction coordinate (
$Q$ ) was defined as$n/9$ (where$n$ is the number of native interactions). The plotted profiles show distinct energy barriers that shift dynamically based on the assigned$E_S$ . - All Interactions vs. Native Only: When non-covalent interactions are permitted between all non-adjacent residues (not just native ones), the unfolding time drastically increases (e.g., mean turns increased from 375 to 743 at E = -1.25 kT).
The simulation outputs the following data for further analysis:
- Total Energy matrices
- Radius of Gyration arrays
- End-to-End distance tracking
- Turn counts for unfolding events
- Microstate frequency distributions
- Animation files (.mp4) generated via
matplotlib.animation